{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "5074e9a3", "metadata": { "tags": [ "hide-output" ] }, "outputs": [], "source": [ "import pymc as pm\n", "import numpy as np\n", "import arviz as az\n", "from itertools import combinations\n", "\n", "%load_ext lab_black" ] }, { "cell_type": "markdown", "id": "d0d4ac84", "metadata": {}, "source": [ "# 7. Coagulation*\n", "\n", "An example of Bayesian ANOVA.\n", "\n", "Adapted from [Unit 7: anovacoagulation.odc](https://raw.githubusercontent.com/areding/6420-pymc/main/original_examples/Codes4Unit7/anovacoagulation.odc)." ] }, { "cell_type": "markdown", "id": "195a19bd", "metadata": {}, "source": [ "Here 24 animals are randomly allocated to 4 different diets, but the numbers allocated to different diets are not the same. The coagulation time for blood is measured for each animal. Are the diet-based differences significant?\n", "\n", "\n", "\n", "[Box, Hunter, Hunter; Statistics for Experimenters, p. 166](https://pages.stat.wisc.edu/~yxu/Teaching/16%20spring%20Stat602/%5bGeorge_E._P._Box,_J._Stuart_Hunter,_William_G._Hu(BookZZ.org).pdf)" ] }, { "cell_type": "code", "execution_count": 2, "id": "3ace28c7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{1: [62, 60, 63, 59],\n", " 2: [63, 67, 71, 64, 65, 66],\n", " 3: [68, 66, 71, 67, 68, 68],\n", " 4: [56, 62, 60, 61, 63, 64, 63, 59]}" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# cut and pasted data from .odc file\n", "# fmt: off\n", "times = (62, 60, 63, 59, 63, 67, 71, 64, 65, 66, 68, 66, 71, 67, 68, 68, 56, 62,\n", " 60, 61, 63, 64, 63, 59)\n", "diets = (1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4)\n", "# fmt: on\n", "\n", "# create dictionary where each key is a diet and values are lists of times\n", "data = {}\n", "for key, val in zip(diets, times):\n", " data.setdefault(key, []).append(val)\n", "data" ] }, { "cell_type": "markdown", "id": "fbfa5f14", "metadata": {}, "source": [ "## Simple method\n", "\n", "No loops! If you're using this style, 4 treatments is probably the max before it starts to get too annoying to type out." ] }, { "cell_type": "code", "execution_count": 3, "id": "47256140", "metadata": { "tags": [ "hide-output" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [μ0, τ, α4, α3, α2]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [24000/24000 00:03<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 4 seconds.\n" ] } ], "source": [ "with pm.Model() as m:\n", " mu0 = pm.Normal(\"μ0\", mu=0, tau=0.0001)\n", " tau = pm.Gamma(\"τ\", 0.001, 0.001)\n", "\n", " alpha4 = pm.Normal(\"α4\", mu=0, tau=0.0001)\n", " alpha3 = pm.Normal(\"α3\", mu=0, tau=0.0001)\n", " alpha2 = pm.Normal(\"α2\", mu=0, tau=0.0001)\n", " # sum-to-zero constraint\n", " alpha1 = pm.Deterministic(\"α1\", -(alpha2 + alpha3 + alpha4))\n", "\n", " mu_1 = mu0 + alpha1\n", " mu_2 = mu0 + alpha2\n", " mu_3 = mu0 + alpha3\n", " mu_4 = mu0 + alpha4\n", "\n", " pm.Normal(\"lik1\", mu=mu_1, tau=tau, observed=data[1])\n", " pm.Normal(\"lik2\", mu=mu_2, tau=tau, observed=data[2])\n", " pm.Normal(\"lik3\", mu=mu_3, tau=tau, observed=data[3])\n", " pm.Normal(\"lik4\", mu=mu_4, tau=tau, observed=data[4])\n", "\n", " onetwo = pm.Deterministic(\"α1-α2\", alpha1 - alpha2)\n", " onethree = pm.Deterministic(\"α1-α3\", alpha1 - alpha3)\n", " onefour = pm.Deterministic(\"α1-α4\", alpha1 - alpha4)\n", " twothree = pm.Deterministic(\"α2-α3\", alpha2 - alpha3)\n", " twofour = pm.Deterministic(\"α2-α4\", alpha2 - alpha4)\n", " threefour = pm.Deterministic(\"α3-α4\", alpha3 - alpha4)\n", "\n", " trace = pm.sample(5000)" ] }, { "cell_type": "code", "execution_count": 4, "id": "a1f8ef11", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_3%hdi_97%
α4-2.9950.817-4.569-1.499
α33.9940.8772.3785.707
α22.0030.8850.3283.661
α1-3.0021.008-4.943-1.137
α1-α2-5.0051.589-8.087-2.071
α1-α3-6.9971.582-9.985-4.006
α1-α4-0.0071.507-2.8152.912
α2-α3-1.9911.419-4.6650.697
α2-α44.9981.3552.4327.570
α3-α46.9891.3424.4889.550
\n", "
" ], "text/plain": [ " mean sd hdi_3% hdi_97%\n", "α4 -2.995 0.817 -4.569 -1.499\n", "α3 3.994 0.877 2.378 5.707\n", "α2 2.003 0.885 0.328 3.661\n", "α1 -3.002 1.008 -4.943 -1.137\n", "α1-α2 -5.005 1.589 -8.087 -2.071\n", "α1-α3 -6.997 1.582 -9.985 -4.006\n", "α1-α4 -0.007 1.507 -2.815 2.912\n", "α2-α3 -1.991 1.419 -4.665 0.697\n", "α2-α4 4.998 1.355 2.432 7.570\n", "α3-α4 6.989 1.342 4.488 9.550" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "az.summary(trace, var_names=[\"α\"], filter_vars=\"like\", kind=\"stats\")" ] }, { "cell_type": "code", "execution_count": 5, "id": "a8f1e0d4", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([], dtype=object)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAKECAYAAACThZlNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuOElEQVR4nO3deXRUdZr/8U8VIQmEgAQSAsgS1oEEQTAwIRBwAe2mZzrIKBDQoCyiMtjYjo3s2hBwAZvWRhAUR9xaDEg3SIMsBmW6WToI2EMCaIIeEWITliCEbPf3B7/U+Mgila1S8f06h3PIrapvPVVoeHPvrRuX4ziOAAAA/j+3rwcAAADVC3EAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAPCBY8eOadKkSWrfvr2Cg4PVuHFj3XHHHdqwYYNX6yxbtkwul0sul0tjxowp8zwHDhzQiBEj1LRpUwUHB6tt27Z67LHHdOrUqcvev7i4WDNmzFCLFi0UFBSkG264QatWrbri+nv37lVAQID+67/+q0zzffTRR57XeTXZ2dme+2VnZ5vbRo0a5bmt9FfdunXVtGlT/eu//qsmTJigzZs362pXlJ81a5ZcLpf69+9fptcB+A0HQJXat2+f06RJE0eSExQU5PTo0cNp166dI8mR5MydO/ea1snJyXHCwsI8jxs9enSZ5tmyZYtTp04dR5ITHh7udO/e3albt64jyWnTpo1z7NixSx7z+OOPO5Kc0NBQp0uXLk5AQIDjcrmcNWvWXPY5+vTp4zRt2tQ5c+ZMmWbcunWr53VeTVZWlud+WVlZ5rbk5GRHkhMREeHEx8c78fHxTq9evZyOHTs6tWvX9jyua9euzr59+y67/syZMx1JTr9+/cr0OgB/QRwAVaiwsNDp0KGDI8np37+/k5OT47lt8+bNTmhoqONyuZy0tLQfXWvEiBGO2+12Bg0aVOY4OHPmjBMeHu5IciZOnOgUFBQ4juM4//znP534+HhHkjNo0CDzmG+//dYJDg52WrVq5QmHtLQ0x+12O127dr3kOV5//XVHkrNixQqv5ytVkXGQnJx8yePOnTvnpKamOl26dHEkOXXr1nXS09MvuR9xgJ8KDisAVWjdunU6ePCggoKC9Nprryk8PNxz2y233KKpU6fKcRw9+eSTV11n06ZNevPNN/XAAw/opptuKvM8ixcv1rfffqtOnTppwYIFql27tiSpUaNGeuuttxQQEKB169YpPT3d85j9+/crPz9f9913n5o0aSJJSkhIUJ8+fbR3717l5eV57puXl6ff/OY36tOnj0aOHFnmOStbnTp1dOedd2rHjh269dZbde7cOd19990qLi729WiATxAHQBXavn27JCk2NlatWrW65PYhQ4ZIuniMPScn57Jr5Ofn68EHH1RERIRSUlLKNU/peQKjRo1SrVq1zG0tW7bUbbfdJkl67733PNtL5yoNg1JNmzaVJJ05c8azbdasWcrJydELL7xQrjmrSp06dfTGG28oKChIhw8f1sqVK309EuATxAFQhU6ePClJat68+WVvL91eUlKiXbt2XfY+s2fP1uHDh/Xss8/quuuuK/MsRUVF+vvf/y5Jio+Pv+x9Srfv2LHDs61ly5aSpIMHD5r7ZmZmKiAgQI0aNZJ08STHF154QQ888IC6detW5jmrWmRkpBITEyVd3NMD/BQRB0AVatCggSTp66+/vuzt39+emZl5ye0HDhzQs88+q759++ree+8t1yzZ2dkqLCyUJLVp0+ay9yndfujQIc+2rl27KiIiQq+88oo2bdqkvLw8LVy4UJ9++qkSEhIUHBwsSfrP//xPNWjQQLNnzy7XnL7Qp08fSbpioAE1XYCvBwB+SmJjYyVJu3fv1ldffaUWLVqY27//ccDSvQylHMfRAw88oJKSEi1atKjcs3x//YYNG172PqXbv3/funXrau7cuRo9erQGDBjg2V6vXj3Nnz9fkrRy5Upt3rxZL7/8smeNwsJC/fOf/1SjRo0UGBhY5rl/7OOMFaH0z+VKh3aAmo44AKrQL3/5SzVr1kxHjx5VUlKS3n33Xc+x+nXr1mnOnDme+54/f9489pVXXtHHH3+sxx57TDExMeWeJT8/3/P7K/1lHRQUdNlZ7r//fjVr1kzLly/Xt99+qw4dOmjSpEnq2LGjzp07p8cee0w33XSTRo8eLcdxNG3aNC1cuFDfffedQkJCNHHiRM2ZM6dMf9Ff6RCIJF24cEG7d+/2es0fCgkJkSRzciXwU0IcAFUoODhYf/zjH/Xzn/9cn3zyiVq2bKmOHTvq5MmTOnr0qFq2bKlu3bpp27Ztqlevnudx3377rX7zm9/o+uuv18yZMytsllIFBQXm61IXLlyQdPFEvR+64447dMcdd1yyfc6cOfrqq6/07rvvyu12a/bs2UpJSdEvfvEL/cd//IdWrVqluXPnKiQkRFOnTvV67k8++eSKt2VnZysqKsrrNX/o7NmzkqT69euXey3AH3HOAVDF+vTpo/T0dN1///2KjIz0nNg3fvx47d692/PxucjISM9jHn/8ceXm5ur555830VAe3z+U8MNDGD/cfqXDDj/0+eefa/78+Ro1apR69eqlwsJCzZ8/X+3atdOaNWuUnJys1atXq127dpo/f76KiorK/0IqwZdffilJioiI8PEkgG+w5wDwgXbt2umVV165ZHtRUZH27t0rSerRo4dn+549eyRJEyZM0IQJE8xjSv+V+9Zbb2nt2rWSLl6e+ce0bt1atWvXVmFhob744gvP4Y3v++KLLyRJ7du3v5aXpUceeUTBwcGaN2+eJCkjI0OnTp1SUlKS3O6L/xZxu90aOHCgFi1apMzMTEVHR1/T2lWpdO9Ez549fTwJ4BvEAVCNbNiwQWfPnlWzZs3UvXv3S24/fvz4FR97/vz5S84NuJqAgAB1795dO3bs0Pbt2y97LL/0ugy9evX60fXWrl2rdevWaeHChZ5/cZeGS2hoqLlv6ddX+tkNvvTNN9/oT3/6kyRp0KBBPp4G8A0OKwDVREFBgWbMmCFJevDBB81FiT799FM5Fy93fsmv0nMQSk/+c67yg4N+6M4775Qkvfbaa5dcDfDLL7/Upk2bJP3fxZmu5MKFC/rVr36lmJgYPfTQQ57tpWf9f/755+b+pV83btz4mmetCufPn9c999yjCxcuqEOHDj/6uoGaijgAqtgHH3xgLiokSV999ZUSExOVnp6uzp07l/mnF17Oe++9p9atW3s+u/9948ePV+PGjXXgwAE9+uijnusenDhxQklJSSoqKtLPfvYzc4jjcp555hl9/vnnevHFFxUQ8H87JJs3b64WLVroz3/+s/bt2yfp4uWX//znPysyMvKaD1dUtvPnz2v16tXq1auXNm/erJCQEL377ruXXDUS+KngsAJQxTZu3KiFCxeqYcOGat26tfLz85WRkSHHcdS5c2dt3LjR8xHCinD27FkdOXLksrfVr19f77zzjn7xi1/o97//vd5++221bNlSBw4c0Llz59S6dWu9+uqrV13/yy+/1Lx58zRs2DD169fP3OZyuTRr1iyNHj1asbGx6tixow4ePKgLFy5o5syZnvMQqtL69es9oVRcXKyTJ0/qiy++8IRRt27dtGLFigr5uCjgr4gDoIolJibqm2++0c6dO3XgwAEFBQUpNjZWQ4cO1cMPP1yhYXAtbr31Vu3evVuzZ8/Wli1btH//fjVv3lyDBw/WtGnTfvSTCo8++qhcLpeee+65y95+//33Kz8/X88//7wyMjLUqlUr/frXv9b48eMr4+X8qJycHM/FjYKDg9WgQQN1795dN910kwYPHqxbb73VJ3MB1YnL8eYAJQAAqPE45wAAABjEAQAAMIgDAABgEAcAAMAgDgAAgOF3H2V0HIcfowoAQDmEhoZe9Uem+10c5OXlqUGDBr4eAwAAv3X69Omr/khyv7vOAXsOAAAonx/bc+B3cQAAACoXJyQCAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMCo9Dh44YUXdN999+mGG25QQECAXC6XPvroo8p+WgAAUEYBlf0EEydOlCQ1bdpU4eHhOnbsWGU/JarAoeN5yj5xTq0b1VX7JqG+HgcAUIEqPQ7Wrl2rHj16KDIyUuPHj9eSJUsq+yl/8k6cvVBpa588V6Cpqz/Tjqxcz7ZeUWGaMzhGDesGVspzNqoXVCnrAgAuz+s4KC4u1oIFC7Rs2TIdOXJEzZo1U3JysqZNm6bMzExFR0dryZIlGjdunCRp0KBBFT50ZThXUOTrESpMj9mbKm1tt0sKDa6tPyR1V2xUQ+3KOqkpq/dr4PPbVOJUznP+71O3V87CVaRuYKU3OABUKK+/aw0dOlSpqanq3bu3EhMTlZaWplmzZqmwsFABAQFyu91KTEyshFErV+cZG3w9gl8ocaSUwV006IamkqRBNzSVI0cT3tpTac/p73822fP8I5ABoJRXcbBq1SqlpqYqISFBW7duldvtVklJieLi4vTSSy8pPDxc8fHxioiIqKx5UQ3ERjU0X/eMCvPRJACAyuBVHKxcuVKSNGnSJLndFz/o4Ha7NWbMGI0bN065ubl68MEHK37KKuDvu66/r7L/pb0r66Rnz4Ek7fze+QeVoSb92QCAP/AqDg4cOCBJio+PN9vj4uI8vx88eHAFjFX1atJx4b9Pu63S1n7ozXRNX/OZHDnqGRWmnVm5mrHmH+oVFaZFI7pXynPWpD8bAPAHXn3XPXv2rGrVqqXw8HCzvVOnTgoODlZMTIxatmxZoQPCe5V5dv9LI3vokXf2mHMM+rZvrIXDblRYSOV8WgEAULW8ioOQkBAVFxcrPz9fwcHBnu1ZWVnKz89X8+bNK3xAVC9hIYFaMboX1zkAgBrMqyskdurUSZKUnp5utq9YsUKStHfv3goaC9Vd+yahGtC5CWEAADWQV3FQ+hHFWbNmqaCgQJKUm5vrubBRdna2MjMzK3ZCAABQpVyO41zzpWtKSkqUkJCg7du3q2vXrrrlllu0fv16ZWRkKCUlRVOnTlVMTIySkpI0efJkSdK8efOUkZEhSfrrX/+qgwcP6vbbb1dkZKQkacyYMerTp08lvDQAAFAWXsWBJJ0+fVozZ85Uamqqjh8/rrCwME2ZMkUTJ07U0qVLNX36dOXn5+vUqVOSpP79+ystLe2K6y1fvlyjRo0qz2sAAAAVyOs4AAAANVul/8hmAADgX4gDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAIxKjYNDhw4pJSVFCQkJatasmQIDA9WiRQvde++9ysjIqMynBgAAZeRyHMeprMWHDRumP/7xj4qJiVGfPn1Uv3597d+/X+vXr1edOnW0YcMG9e3bt7KeHgAAlEGlxsFrr72mG2+8UV27djXb33nnHQ0fPlydO3fWP/7xj8p6egCocoeO5yn7xDm1blRX7ZuE+nocoEy8joPi4mItWLBAy5Yt05EjR9SsWTMlJydr2rRpyszMVHR0tJYsWaJx48ZddZ2OHTvq4MGD+vbbb9W4ceNyvQgA8NaJsxcqdL2T5wo0dfVn2pGV69nWKypMcwbHqGHdwAp7nkb1gipsLeBKArx9wNChQ5WamqrevXsrMTFRaWlpmjVrlgoLCxUQECC3263ExMQfXad27doXBwjwegQA1dC5giJfj+CVHrM3Veh6bpcUGlxbf0jqrtiohtqVdVJTVu/XwOe3qaQC98/+71O3V9xilaxuIN/f/ZVXew5WrVqlIUOGKCEhQVu3bpXb7VZJSYni4uJ0+PBhhYeHKyIiQtu2bbvqOjt37lSvXr0UGxurnTt3lvtFAPC91pPX+XoEn/tDUncNuqGp5+u1+45qwlt7fDiRb2XPG+TrEVBGXn1aYeXKlZKkSZMmye2++FC3260xY8YoNzdXmZmZGjJkyFXXOH36tJKTk+V2u/XMM8+UcWwAqH5ioxqar3tGhfloEqB8vNrnc+DAAUlSfHy82R4XF+f5/eDBg6/4+Pz8fN15553KyMjQnDlz1L9/f2+eHkA15k+7uyWp84wNFb7mrqyTZs/Bzu+df1BR/O19hn/y6rBCu3btlJ2draIie2yxuLhY9erVU0xMjHbt2nXZx164cEGJiYn6y1/+oieeeEIpKSnlmxwAyqGiT0h86M10Hco5q6d+Ga2eUWHamZWrGWv+ofYR9bRoRPcKex5OSERV8GrPQUhIiIqLi5Wfn6/g4GDP9qysLOXn56t58+aXfVx+fr4SExO1YcMGPf7444QBAJ+r6L9kXxrZQ4+8s8ecY9C3fWMtHHajwkIq7tMKQFXwKg46deqkffv2KT09Xb179/ZsX7FihSRp7969lzzm+2Hw2GOP6emnny7nyABQ/YSFBGrF6F5c5wA1gleHFUovXjRgwACtXbtWgYGBys3NVefOnXX8+HFJUkZGhjp27CjpYhj88pe/1MaNG/Xoo49q/vz5lfMqAABAhfEqDkpKSpSQkKDt27era9euuuWWW7R+/XplZGQoJSVFU6dOVUxMjJKSkjR58mSNGjVK//3f/63IyEg98MADl11z1KhRat26dUW9HgAAUE5eXyHx9OnTmjlzplJTU3X8+HGFhYVpypQpmjhxopYuXarp06crPz9fp06dUv/+/ZWWlnbV9bZu3cqnFgAAqEYq9WcrAAAA/1OpP7IZAAD4H+IAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAARqXGwZEjRzR+/Hj16NFD4eHhCgoKUqtWrTRo0CBt3ry5Mp8aAACUkctxHKeyFt+0aZPuuusuxcXFKSoqSvXr19fXX3+tNWvW6MyZM5ozZ46mTJlSWU8PAPAjh47nKfvEObVuVFftm4T6epyftEqNg4KCAgUEBMjttjsojh49qu7duys3N1c5OTm67rrrKmsEAEAFO3H2QoWud/Jcgaau/kw7snI923pFhWnO4Bg1rBtYYc/TqF5Qha1V03kdB8XFxVqwYIGWLVumI0eOqFmzZkpOTta0adOUmZmp6OhoLVmyROPGjbvqOnfeeadWr16tTz/9VF27di3XiwCAmuJcQZGvR/hRnWdsqND13C4pNLi2UgZ3UWxUQ+3KOqkpq/crL79QJRX4z9f/fer2ilusitQNDPDJ83r9rEOHDlVqaqp69+6txMREpaWladasWSosLPTsJUhMTLzqGidOnNCOHTtUt25dtWnTpqyzA0CNU9F/8fqDEkdKGdxFg25oKkkadENTOXI04a09Ffo8/vjeZs8b5JPn9SoOVq1apdTUVCUkJGjr1q1yu90qKSlRXFycXnrpJYWHhys+Pl4RERHmcdnZ2XrttddUXFyso0eP6k9/+pNOnTqlxYsXKzSU40oA8FMXG9XQfN0zKsxHk0DyMg5WrlwpSZo0aZLnPAK3260xY8Zo3Lhxys3N1YMPPnjJ47Kzs/Xkk096vq5Xr56WL1+ukSNHlmd2AKhx/GHXd2X8C3xX1knPngNJ2vm98w8qij+8t9WFV3Fw4MABSVJ8fLzZHhcX5/n94MGDL3lc//795TiOCgsLlZ2drZdffln33nuvdu7cqd///vdlmRsAaiRfHWP2xt+n3Vah6z30Zrqmr/lMjhz1jArTzqxczVjzD/WKCtOiEd0r7Hn84b2tLrw6IbFdu3bKzs5WUZE9Yaa4uFj16tVTTEyMdu3adU1rPfzww1q0aJE++OAD/exnP/NuagBAjZH7XYEeeWePPj70T8+2vu0ba+GwGxUWUnGfVsC18yqjQkJCVFxcrPz8fAUHB3u2Z2VlKT8/X82bN7/mtQYOHKhFixbpo48+Ig4A4CcsLCRQK0b34joH1YhXV0js1KmTJCk9Pd1sX7FihSRp796917zW0aNHJUkBAezmAQBI7ZuEakDnJoRBNeBVHJR+RHHWrFkqKCiQJOXm5mrJkiWSLp54mJmZ6bn/zp07lZ+ff8k6R44c0dy5cyWJvQYAAFQzXp1zUFJSooSEBG3fvl1du3bVLbfcovXr1ysjI0MpKSmaOnWqYmJilJSUpMmTJysxMVEff/yx+vXrp5YtWyogIECff/65PvjgAxUUFGjSpElasGBBZb4+AADgJa+vkHj69GnNnDlTqampOn78uMLCwjRlyhRNnDhRS5cu1fTp05Wfn69Tp05p7dq1evPNN7Vr1y4dO3ZMBQUFioiIUM+ePTV27Fj2GgAAUA1V6s9WAAAA/qdSf2QzAADwP8QBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAEaVx8FDDz0kl8sll8ulY8eOVfXTAwCAH1GlcbB582YtXrxYISEhVfm0AADACy7HcZyqeKK8vDx16dJFPXr00IkTJ5SWlqZvvvlGkZGRVfH0qGCHjucp+8Q5tW5UV+2bhPp6HABABQrw9gHFxcVasGCBli1bpiNHjqhZs2ZKTk7WtGnTlJmZqejoaC1ZskTjxo0zj/v1r3+tvLw8LVq0SEOHDq2wF4Brd+LshXKvcfJcgaau/kw7snI923pFhWnO4Bg1rBtY7vUb1Qsq9xoAgPLxOg6GDh2q1NRU9e7dW4mJiUpLS9OsWbNUWFiogIAAud1uJSYmmsds3LhRS5cu1euvv64mTZpU1OzVxrmCIl+PcE16zN5U7jXcLik0uLb+kNRdsVENtSvrpKas3q+Bz29TSQXsg/rfp24v/yKVpG6g1/+7AIBf8uqwwqpVqzRkyBAlJCRo69atcrvdKikpUVxcnA4fPqzw8HBFRERo27ZtnsecOXNGMTEx6tKli9atWydJ6t+/f406rNB68jpfj1Cl/pDUXYNuaOr5eu2+o5rw1h4fTlQ1sucN8vUIAFAlvDohceXKlZKkSZMmye2++FC3260xY8YoNzdXmZmZGjJkiHnMr371K50+fVpLliypoJHha7FRDc3XPaPCfDQJAKAyeLWf9MCBA5Kk+Ph4sz0uLs7z+8GDB3t+v379ei1fvlyLFy/W9ddfX545q7XqvCv8+zrP2FAh6+zKOmn2HOz83vkH5eUv7yUA1GRexcHZs2dVq1YthYeHm+2dOnVScHCwYmJi1LJlS0nSuXPnNHbsWN18882XnJxY0/jLsei/T7ut3Gs89Ga6pq/5TI4c9YwK086sXM1Y8w/1igrTohHdy72+v7yXAFCTefWdOCQkRMXFxcrPz1dwcLBne1ZWlvLz89W8eXPPtpycHH399df6+uuvPYcgfqhp04v/+tyzZ4+6detWhvHhjYr4JMBLI3vokXf2mHMM+rZvrIXDblRYSPk/rQAA8D2v4qBTp07at2+f0tPT1bt3b8/2FStWSJL27t3r2RYaGqrRo0dfdp1169bp2LFjSkpKUp06ddSoUaOyzA4fCAsJ1IrRvbjOAQDUYF59WuGdd97R8OHDNWDAAK1du1aBgYHKzc1V586ddfz4cUlSRkaGOnbseNV1atqnFQAAqEm8+rTC3Xffrfj4eH344Yfq2bOnHn30UcXHx+v48eNKSUmRy+XSXXfdpXnz5lXWvAAAoJJ5FQdut1vr1q3TI488ohMnTujFF1/UyZMntXDhQj3xxBNasmSJcnJyiAMAAPxYlf1sBQAA4B+q/Ec2AwCA6o04AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwiAMAAGAQBwAAwCAOAACAQRwAAACDOAAAAAZxAAAADOIAAAAYxAEAADCIAwAAYBAHAADAIA4AAIBBHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMQBAAAwAnw9gLccx1FeXp6vxwAAwG+FhobK5XJd8Xa/i4O8vDw1aNDA12MAAOC3Tp8+rfr161/xdpfjOE4VzlNuP6U9B2fOnFGLFi301VdfXfUPEZfivSs73rvy4f0rO967svP2vatxew5cLtdP7j+a+vXr/+Rec0XhvSs73rvy4f0rO967squo944TEgEAgEEcAAAAgzioxoKCgjRz5kwFBQX5ehS/w3tXdrx35cP7V3a8d2VX0e+d352QCAAAKhd7DgAAgEEcAAAAgzgAAAAGcQAAAAziAAAAGMSBn/rb3/6mWrVqyeVyad68eb4ep1r77rvv9MYbb+juu+9Whw4dVKdOHV133XXq16+f3n77bV+PV23s2rVLP//5z9WwYUOFhISoZ8+eeuutt3w9VrX39ddf63e/+50GDhyoli1bKjAwUJGRkRoyZIh27Njh6/H8zjPPPCOXyyWXy6W//e1vvh7Hb6xevVoDBgxQo0aNVKdOHUVFRWn48OH66quvyrSe310+GdL58+c1atQo1alTR999952vx6n2Pv74Y91zzz1q1KiRbr31Vg0ZMkQ5OTlatWqVkpKS9D//8z964YUXfD2mT3300Ue6/fbbFRgYqGHDhqlBgwZatWqVRowYoezsbE2ZMsXXI1ZbL7zwgp5++mm1bdtWAwYMUEREhA4dOqT3339f77//vt5++23dfffdvh7TLxw4cEAzZsxQSEgI39uukeM4Gj9+vF5++WW1bdtWw4YNU2hoqI4ePaq0tDQdOXJELVq0KNPC8DOTJk1y6tev7/z2t791JDlz58719UjV2qeffuq8+eabTkFBgdl+7Ngxp1WrVo4kZ+fOnT6azvcKCwudtm3bOkFBQU56erpn+5kzZ5zo6GgnICDAOXjwoA8nrN5SU1Odbdu2XbJ927ZtTu3atZ2wsDAnPz/fB5P5l6KiIic2Ntbp2bOnM3LkSEeS89e//tXXY1V7CxcudCQ5Dz/8sFNUVHTJ7YWFhWVal8MKfmb79u1auHChnnvuOV1//fW+HscvdO3aVUlJSapdu7bZ3qRJEz3wwAOSpLS0NF+MVi1s2bJFn3/+uZKSknTjjTd6toeGhmr69OkqKirS8uXLfThh9XbnnXeqb9++l2zv27evbr75ZuXm5mr//v0+mMy/PP3009q7d69effVV1apVy9fj+IXz58/rySefVJs2bfS73/3usu9bQEDZDhAQB37k3LlzGjVqlPr376+xY8f6epwaoTQYyvo/UE3w0UcfSZIGDhx4yW2l237K8VQe/Pd1bT777DM9+eSTmjZtmqKjo309jt/48MMPlZubq8TERBUXF2vVqlWaN2+eFi9erMOHD5drbf6L9SOTJ0/WN998o40bN/p6lBqhuLhYr7/+ulwul2677TZfj+Mzhw4dkiS1b9/+ktsaNmyoxo0be+6Da/fll19q06ZNioyMVJcuXXw9TrVVVFSkUaNGqVOnTpo8ebKvx/Eru3fvlnQxPrt27arMzEzPbW63W5MmTdJzzz1XprXZc+An0tLS9OKLLyolJUVRUVG+HqdGmD59uvbv36/77rtPMTExvh7HZ06fPi1JatCgwWVvr1+/vuc+uDaFhYW65557dOHCBT3zzDPsJr+KlJQUz+GEHx76w9Xl5ORIkubPn6/69etr586dysvL07Zt29ShQwfNnz9fL730UpnWJg6qUOPGjT0f0bmWX6W7e7/77jvdf//9iouL04QJE3z7InyorO/f5bz88suaO3eubrzxRi1cuLDqXgRqvJKSEt1///3atm2bxo4dq3vuucfXI1Vbe/fu1ezZs/XYY4+pe/fuvh7H75SUlEiSAgMD9f777ys2Nlb16tVT37599d5778ntdmv+/PllWpvDClVo+PDhysvLu+b7R0ZGSpKmTp2qo0eP6oMPPpDb/dPtubK+fz+0fPlyjR8/Xl26dNGHH36oevXqVdSIfql0j8GV9g6cOXPminsVYDmOo7Fjx+qNN97QyJEjtXjxYl+PVK0lJyerbdu2mjVrlq9H8Uul/1/edNNNatasmbktOjpabdq00eHDh3Xq1Cldd9113i1ezk9RoAr069fPkfSjvx555BFfj1rtvfLKK47b7Xaio6OdnJwcX49TLTzxxBOOJOftt9++5Lbc3FxHktO7d28fTOZfiouLnfvuu8+R5AwfPvyyHyuDdS3f1yQ5q1ev9vWo1dLSpUsdSc6//du/Xfb2m266yZHkHD161Ou12XPgBwYNGqR27dpdsv3QoUPatm2bYmNjdcMNNyguLs4H0/mPV199VWPGjFGnTp20ZcsWhYeH+3qkaqFfv36aO3euNm7cqGHDhpnbSk9+7devny9G8xslJSUaM2aMli9frqFDh2rFihWcZ3ANRo8efdnt27Zt06FDh/Tv//7vCg8PV+vWrat2MD9x8803S7p48agfKiws1OHDhxUSElK273XlLRf4zvLly7kI0jVatmyZ43K5nE6dOjnHjh3z9TjVSmFhodOmTRsnKCjI2bNnj2f79y+ClJmZ6bsBq7ni4mJn1KhRjiTnrrvuKvNFZ/B/kpOTuQjSNRo4cKAjyVm6dKnZ/tRTTzmSnJEjR5ZpXfYcoMbbsmWLxo4dK8dxlJCQcNmzd7t166bExMSqH64aCAgI0LJly3T77berb9++Gj58uOrXr69Vq1YpKytLs2fPVocOHXw9ZrX11FNP6bXXXlO9evXUoUMHzZ49+5L7JCYmqlu3blU/HGq8RYsWqXfv3ho7dqzef/99/cu//Iv27NmjLVu2qFWrVnr22WfLtC5xgBrvyy+/lOM4kqQlS5Zc9j7Jyck/2TiQLu6e/OSTTzRz5ky9++67KigoUHR0tH77299qxIgRvh6vWsvOzpYknT17VnPmzLnsfVq3bk0coFK0bdtWu3fv1owZM/SXv/xFGzduVGRkpB5++GHNmDFDERERZVrX5ZR+1wQAABDXOQAAAD9AHAAAAIM4AAAABnEAAAAM4gAAABjEAQAAMIgDAABgEAcAAMAgDgAAgEEcAAAAgzgAAAAGcQAAAIz/B9naiJZUHFcLAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "az.plot_forest(trace, var_names=[\"α1\", \"α2\", \"α3\", \"α4\"], combined=True)" ] }, { "cell_type": "code", "execution_count": 11, "id": "675ce0bb", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([], dtype=object)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAASBCAYAAAAjXAT/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABfe0lEQVR4nO3de1iVdb7//9finAgkioiaiXn4WuQZFTAlR8W2mmhNmody1ElznEx/NTZZhu20g7vttmZ0cpdSFtp4QsuUJvMwpYkpaqYkOpClk45nCEGEz+8PN2skQFkKLj76fFwX1yX3ab2Xayae3ve9Fg5jjBEAAIDFPNw9AAAAwLUiaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAV/TTTz9p4sSJatasmfz8/FSnTh317t1bKSkpLh3n7bfflsPhkMPh0OjRo696nn379mno0KEKCwuTn5+f7rjjDj311FM6ffp0mdsXFhZq6tSpuu222+Tr66tWrVpp+fLl5R5/165d8vLy0tNPP31V823YsMH5PC8nKyvLuV1WVlaJdSNGjHCuK/6qUaOGwsLC1LlzZ40fP17r1q3T5X57TUJCghwOh2JjY6/qeQBWMQBwGbt37zahoaFGkvH19TXt27c3TZs2NZKMJPPyyy9X6DjHjh0zwcHBzv1GjRp1VfN8/vnn5pZbbjGSTEhIiGnXrp2pUaOGkWSaNGlifvrpp1L7/OEPfzCSTEBAgLn77ruNl5eXcTgcZuXKlWU+RpcuXUxYWJg5e/bsVc24fv165/O8nMzMTOd2mZmZJdY9+uijRpKpW7euiYmJMTExMaZTp06mRYsWxtvb27lf69atze7du8s8/gsvvGAkmW7dul3V8wBswhkaAOW6cOGCHnzwQR09elSxsbH64Ycf9PXXXysjI0Pr1q1TQECAnn32WW3atOmKx5o4caJOnz6tPn36XPU82dnZGjRokM6dO6cnnnhChw8f1vbt23Xo0CHFxMToH//4h0aNGlVin+PHj+uNN97Q7bffroyMDO3evVvr1q2Tw+HQ1KlTSz3GwoUL9cUXX+i1115TQEDAVc9aWe677z598cUX+uKLL/TVV18pPT1dZ86c0bJly3T33Xdr165d6ty5s9LS0tw9KuBWBA2Acq1evVr79++Xr6+vEhMTFRIS4lzXvXt3TZkyRcYYTZs27bLH+eyzz/TBBx9ozJgx6tChw1XP85e//EX/+te/1LJlS/33f/+3vL29JUm1a9dWUlKSvLy8tHr1au3YscO5zzfffKO8vDz95je/UWhoqCSpa9eu6tKli3bt2qXs7GznttnZ2Zo8ebK6dOmiYcOGXfWcVe2WW27RwIEDtXXrVv3qV79Sbm6uHnroIRUWFrp7NMBtCBoA5fryyy8lSZGRkbr99ttLrX/ggQckXbxn5NixY2UeIy8vT48//rjq1q2rGTNmXNM8xfe9jBgxQp6eniXWNWrUSD169JAkLV261Lm8eK7imCkWFhYmSTp79qxzWUJCgo4dO6Y333zzmua8Xm655Ra9//778vX11YEDB7RkyRJ3jwS4DUEDoFynTp2SJDVo0KDM9cXLi4qKtG3btjK3eemll3TgwAHNnDlTt95661XPcuHCBW3fvl2SFBMTU+Y2xcu3bt3qXNaoUSNJ0v79+0ts+91338nLy0u1a9eWdPFG4zfffFNjxoxRmzZtrnrO661evXqKj4+XdPGMGnCzImgAlCsoKEiSdPjw4TLXX7r8u+++K7V+3759mjlzpu655x498sgj1zRLVlaWCgoKJElNmjQpc5vi5RkZGc5lrVu3Vt26dfXOO+/os88+U3Z2tmbPnq2dO3eqa9eu8vPzkyT9/ve/V1BQkF566aVrmtMdunTpIknlRiVwM/By9wAAqq/IyEhJ0tdff60ffvhBt912W4n1l771ufhsTjFjjMaMGaOioiLNmTPnmme59Pi1atUqc5vi5ZduW6NGDb388ssaNWqUevbs6Vxes2ZNvf7665KkJUuWaN26dZo3b57zGAUFBTp+/Lhq164tHx+fq577Sm/drgzFr0t5l/2AmwFBA6Bc/fv3V/369XXkyBENGTJEf/3rX533nqxevVrTp093bnvu3LkS+77zzjv6+9//rqeeekoRERHXPEteXp7zz+UFhq+vb5mzjBw5UvXr19eCBQv0r3/9S82bN9fEiRPVokUL5ebm6qmnnlKHDh00atQoGWP03HPPafbs2fr555/l7++vJ554QtOnT7+qOCnv8pgk5efn6+uvv3b5mL/k7+8vSSVucAZuNgQNgHL5+fnpww8/1H/8x3/oiy++UKNGjdSiRQudOnVKR44cUaNGjdSmTRtt2rRJNWvWdO73r3/9S5MnT1bDhg31wgsvVNosxc6fP1/i+2L5+fmSLt4s+0u9e/dW7969Sy2fPn26fvjhB/31r3+Vh4eHXnrpJc2YMUN9+/bVgw8+qOXLl+vll1+Wv7+/pkyZ4vLcX3zxRbnrsrKyFB4e7vIxfyknJ0eSFBgYeM3HAmzFPTQALqtLly7asWOHRo4cqXr16jlvrh07dqy+/vpr51uF69Wr59znD3/4g06ePKlZs2aVCJ1rcellpl9e3vrl8vIuSf3SwYMH9frrr2vEiBHq1KmTCgoK9Prrr6tp06ZauXKlHn30Ua1YsUJNmzbV66+/rgsXLlz7E6kChw4dkiTVrVvXzZMA7sMZGgBX1LRpU73zzjulll+4cEG7du2SJLVv3965vPhD3saPH6/x48eX2Kf4bEJSUpI+/vhjSRd/tcKVNG7cWN7e3iooKNA//vEP56WvS/3jH/+QJDVr1qwiT0sTJkyQn5+fXnnlFUlSenq6Tp8+rSFDhsjD4+K/9zw8PNSrVy/NmTNH3333ne66664KHft6Kj4L1LFjRzdPArgPQQPgqqWkpCgnJ0f169dXu3btSq0/evRoufueO3eu1L0ul+Pl5aV27dpp69at+vLLL8u8N6X4c3M6dep0xeN9/PHHWr16tWbPnu08s1EcW7/8hODi78v7XVHu9M9//lOrVq2SpGv6FGbAdlxyAnBVzp8/7/zVAY8//niJD7rbuXOnjDFlfhXfU1N8A665zC9X/KWBAwdKkhITE0t9Ku6hQ4f02WefSfr3B/6VJz8/X08++aQiIiI0btw45/LidwsdPHiwxPbF39epU6fCs14P586d0/Dhw5Wfn6/mzZtf8XkDNzKCBsBlffLJJyU+qE6SfvjhB8XHx2vHjh268847r/q3Updl6dKlaty4sfOzVS41duxY1alTR/v27dOkSZOcn0tz4sQJDRkyRBcuXNB9991X4vJXWV577TUdPHhQf/rTn+Tl9e8T1Q0aNNBtt92mjz76SLt375Z08VcnfPTRR6pXr16FL2VVtXPnzmnFihXq1KmT1q1bJ39/f/31r38t9enJwM2ES04ALuvTTz/V7NmzVatWLTVu3Fh5eXlKT0+XMUZ33nmnPv30U+fbpStDTk6Ovv/++zLXBQYGavHixerbt6/eeOMNLVq0SI0aNdK+ffuUm5urxo0ba/78+Zc9/qFDh/TKK69o8ODB6tatW4l1DodDCQkJGjVqlCIjI9WiRQvt379f+fn5euGFF5z31VxPa9asccZdYWGhTp06pX/84x/OmGvTpo0WLlxYKW+NB2xG0AC4rPj4eP3zn/9Uamqq9u3bJ19fX0VGRmrQoEH63e9+V6kxUxG/+tWv9PXXX+ull17S559/rm+++UYNGjTQgAED9Nxzz13xHU6TJk2Sw+HQf/3Xf5W5fuTIkcrLy9OsWbOUnp6u22+/Xf/f//f/aezYsVXxdK7o2LFjzg/M8/PzU1BQkNq1a6cOHTpowIAB+tWvfuWWuYDqxmFcuYANAABQDXEPDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsx+fQXAfGGGVnZ7t7DAAArBUQECCHw1HueoLmOsjOzlZQUJC7xwAAwFpnzpxRYGBguev5YL3rgDM0AABcmyudoSFoAACA9bgpGAAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QiaMnz00Uf6/e9/r5iYGPn7+8vhcCghIcHdYwEAgHJ4uXuA6uj111/Xxo0bFRgYqPr16+vAgQPuHgkAAFwGZ2jK8J//+Z/av3+/Tp8+rf/8z/909zgAgJtExtFs/W3vUWUczXb3KNaptkGzYcMG56WezZs3695771VAQIBCQkI0btw4nTt3TpK0du1a56Wh0NBQTZ48WYWFhc7jJCQkyOFwaMOGDaUeIzExUQ6HQ4mJiSWW33PPPWrWrJkcDkdVPkUAAHQiJ18HjmVr0Ftb1HPWJv32va/Vc9YmDXpriw4cy9aJnHx3j2iFan/JaevWrXr11VcVFxenMWPGaP369Zo7d67Onj2r/v3769FHH9X999+vTp06afXq1XrttdcUGBioKVOmuHt0AACUe/7CZde3f+kzeTikAD9v/XlIO0WG19K2zFN6dsU36jVrk4qMtPfFuMseo4ZPtf9xXuWq/d/A2rVrlZycrP79+0uSCgoK1KFDByUlJSklJUUbN25UZGSkJGnatGlq2rSpZs2apcmTJ8vLq9o/PQDADe7OqSlX3KbISDMG3K0+rcIkSX1ahcnIaHxSWoWOkfVKn2sf1HLV9pJTsdjYWGfMSJK3t7cefPBBGWPUr18/Z8xIUkBAgPr27asTJ07oxx9/dMe4AABclcjwWiW+7xge7KZJ7FTtT2G0bdu21LKwsIsF26ZNm3LXHT58WI0bN67K0QAAuKIrXS4qPvuyLfOU8wyNJKVmnqzwMWBB0AQGBpZaVnwp6XLrCgoKqnYwAAAq4Er3t2x/rofGfbBDz6/cIyOjjuHBSs08qakrv1Wn8GDNGdqOe2Qq4Ib/G/LwuHhV7cKF0jdlnTlz5nqPAwBACbVr+mrusPaasDjNec+MJN3TrI5mD26rYH8fN05njxs+aGrVunhN8vDhw6XWpaWllVoGAMD1Fuzvo4WjOinjaLayTuSqce0aahYa4O6xrFLtbwq+Vh06dJAkvffeeyoqKnIu37Jliz744AN3jQUAQCnNQgPU885QYuYq3PBnaDp37qyoqCh9/vnnioqKUteuXfX9999r1apV6tevn1asWFFqn+TkZCUnJ0uSMjMzncuysrIkSV26dNHo0aOv11MAAABXcMMHjcPh0KpVqzRp0iStXr1a33zzjVq3bq1Vq1bpyJEjZQbNzp079e6775ZYtmvXLu3atcv5PUEDAED14TDGGHcPAQAAcC1u+HtoAADAjY+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9guYXTpw4oXnz5un+++9XkyZN5Ovrqzp16ui+++5TSkqKu8cDAABlcBhjjLuHqE7+8pe/6PHHH1eDBg3UvXt3NWjQQD/++KOWLVumc+fOaebMmXrqqafcPSYAALgEQfMLn3/+uc6dO6f77rtPHh7/PoH13XffqVOnTsrNzVVWVpbq16/vxikBAMClqu0lpw0bNsjhcCghIUGbN2/Wvffeq4CAAIWEhGjcuHE6d+6cJGnt2rWKiYmRv7+/QkNDNXnyZBUWFjqPk5CQIIfDoQ0bNpR6jMTERDkcDiUmJjqXde/eXX369CkRM5LUokULDRo0SAUFBdq8eXOVPGcAAK6XjKPZ+tveo8o4mu3uUSqFl7sHuJKtW7fq1VdfVVxcnMaMGaP169dr7ty5Onv2rPr3769HH31U999/vzp16qTVq1frtddeU2BgoKZMmVLps3h7e0uSvLyq/V8bAAClnMjJ16nc85qyYo+2Zp50Lu8UHqy5w9or2N/HjdNdm2r/k3nt2rVKTk5W//79JUkFBQXq0KGDkpKSlJKSoo0bNyoyMlKSNG3aNDVt2lSzZs3S5MmTKzU8srOztXTpUvn5+emee+6ptOMCAFDZcs9fKHN5+5c+k4dDCvDz1p+HtFNkeC1tyzylZ1d8o/FJO/T2ox3K3K+GT7XPheofNLGxsc6YkS6eJXnwwQe1e/du9evXzxkzkhQQEKC+fftq/vz5+vHHH9W4ceNKm2Ps2LE6evSoXnzxRdWuXbvSjgsAQGW7c2r578otMtKMAXerT6swSVKfVmEyMhqflFbuflmv9KmSOStTtb2Hpljbtm1LLQsLu/gitGnTptx1hw8frrQZnn32WSUlJal379569tlnK+24AAC4Q2R4rRLfdwwPdtMklafan6EJDAwstaz4UtLl1hUUFFTK40+bNk0vv/yyunfvruXLl8vT07NSjgsAQFXZ+2JcmcuLz8BsyzzlPEMjSan/dz/NqvExalq3ZtUPWAWqfdBcq+J3K124UPp64pkzZy6777Rp05SQkKDY2Fh99NFHuuWWW6pkRgAAKlN597xsf66Hxn2wQ8+v3CMjo47hwUrNPKmpK7/VPc3qqFXDW6/voJXohg+aWrUunlYr6xJUWlpaufslJCRo2rRp6tatm1avXq0aNWpU2YwAAFwPtWv6au6w9pqwOE3jk/79M/CeZnU0e3DpWzxscsMHTYcOF+/Yfu+99zR8+HDnGZstW7bogw8+KHOfF154QS+++KLuueceYgYAcEMJ9vfRwlGdlHE0W1knctW4dg01Cw1w91jX7IYPms6dOysqKkqff/65oqKi1LVrV33//fdatWqV+vXrpxUrVpTYPjExUS+++KK8vLzUsWNHzZw5s9QxY2NjFRsbe52eAQAAla9ZaMANETLFbvigcTgcWrVqlSZNmqTVq1frm2++UevWrbVq1SodOXKkVNBkZWVJunjPzeuvv17ucQkaAACqD36XEwAAsF61/xwaAACAKyFoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaH7BGKOnn35asbGxql+/vvz8/BQaGqro6Gi98847KigocPeIAADgFxzGGOPuIaqTCxcuqGbNmurQoYNatmypkJAQnTp1SmvXrlVWVpbi4uL0ySefyMODFgQAoLogaMqQl5cnPz+/EssuXLigXr16af369fr444/Vp08fN00HAAB+qdqeZtiwYYMcDocSEhK0efNm3XvvvQoICFBISIjGjRunc+fOSZLWrl2rmJgY+fv7KzQ0VJMnT1ZhYaHzOAkJCXI4HNqwYUOpx0hMTJTD4VBiYmKJ5b+MGUny8vJSfHy8JOnAgQOV9jwBAKjOMo5m6297jyrjaLa7R7ksL3cPcCVbt27Vq6++qri4OI0ZM0br16/X3LlzdfbsWfXv31+PPvqo7r//fnXq1EmrV6/Wa6+9psDAQE2ZMqVS5ygqKtLatWslSREREZV6bAAAqpsDx7I1ZcUebc086VzWKTxY0wdEqGndADdOVrZqHzRr165VcnKy+vfvL0kqKChQhw4dlJSUpJSUFG3cuFGRkZGSpGnTpqlp06aaNWuWJk+eLC+va3t6CQkJkqTjx49r3bp1Sk9P14gRI/SrX/3qmo4LAEB1kHv+Qrnres3apAA/b/15SDtFhtfStsxTenbFN+o1a5P2TIsrc58aPu7LimofNLGxsc6YkSRvb289+OCD2r17t/r16+eMGUkKCAhQ3759NX/+fP34449q3LjxNT32tGnTnH92OBx66qmn9PLLL1/TMQEAqC7unJpy2fUzBtytPq3CJEl9WoXJyGh8Ulq5+2W94r77S6vtPTTF2rZtW2pZWNjFv9w2bdqUu+7w4cPX/NjGGBUWFuqHH37QnDlz9Pbbbys2NlZnz5695mMDAFDdRYbXKvF9x/BgN01yZdX+DE1gYGCpZcWXki63rrI+L8bDw0MNGzbU2LFjVbt2bT300EOaPn26Xn311Uo5PgAA7rL3xbIvHR04lqP7//SltmWecp6hkaTU/7ufZtX4GDWtW/O6zFhR1T5orlXx58VcuFD6OuGZM2dcOlavXr0kqcx3TAEAYJvy7nlp1fBWdQoP1vMr98jIqGN4sFIzT2rqym/VKTxYrRreen0HrYAbPmhq1bp4uqysS1BpaWkuHevIkSOSdM03GwMAUN3NHdZeExanaXzSv39W3tOsjmYPLn0rSHVww/9k7tChgyTpvffe0/Dhw51nbLZs2aIPPvig1Pbp6ekKDg5W3bp1SyzPzc3VpEmTJEn33XdfFU8NAIB7Bfv7aOGoTso4mq2sE7lqXLuGmoVWv7drF7vhg6Zz586KiorS559/rqioKHXt2lXff/+9Vq1apX79+mnFihUltl+7dq0mT56s2NhYNWnSREFBQTp8+LDWrFmjEydOKCYmxhk2AADc6JqFBlTrkCl2wweNw+HQqlWrNGnSJK1evVrffPONWrdurVWrVunIkSOlgqZHjx4aNWqUvvjiC23btk3Z2dkKCgpSRESEBg8erNGjR3PJCQCAaobf5QQAAKxX7T+HBgAA4EoIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CJoKWLJkiRwOhxwOhxYvXuzucQAAwC8QNFdw7NgxjRs3Tv7+/u4eBQAAlIOguYLHH39c/v7+Gjt2rLtHAQAA5ai2QbNhwwY5HA4lJCRo8+bNuvfeexUQEKCQkBCNGzdO586dkyStXbtWMTEx8vf3V2hoqCZPnqzCwkLncRISEuRwOLRhw4ZSj5GYmCiHw6HExMQyZ0hKStLy5cs1b9481axZsyqeJgAAThlHs/W3vUeVcTTb3aNYp9oGTbGtW7fqV7/6lYKCgjRmzBg1atRIc+fO1W9/+1stWbJEAwcO1G233aYxY8bo1ltv1WuvvaZXXnnlmh/3p59+0u9//3uNHDlSvXr1qoRnAgBA2Q4cy9agt7ao56xN+u17X6vnrE0a9NYWnfz5vLtHs4aXuwe4krVr1yo5OVn9+/eXJBUUFKhDhw5KSkpSSkqKNm7cqMjISEnStGnT1LRpU82aNUuTJ0+Wl9fVP70xY8bIz89Pr7/+eqU8DwAAcs9fKHN5r1mbFODnrT8PaafI8FralnlKz674RuOTdujtRzuU2r6GT7X/8X3dVfu/kdjYWGfMSJK3t7cefPBB7d69W/369XPGjCQFBASob9++mj9/vn788Uc1btz4qh7zvffe06pVq7Ry5Urdeuut1/gMAAC46M6pKeWumzHgbvVpFSZJ6tMqTEZG45PSytwn65U+VTajrar9Jae2bduWWhYWdvEFb9OmTbnrDh8+fFWPd+TIET355JMaPHiw7r///qs6BgAArooMr1Xi+47hwW6axE7V/gxNYGBgqWXFl5Iut66goOCqHm/cuHHy9PTUm2++eVX7AwBQnr0vxpVaduBYju7/05falnnKeYZGklIzT0qSVo2PUdO6vDHlSqp90FwrD4+LJ6EuXCh93fLMmTOllu3cuVPHjx9XSEhImcd7+OGH9fDDD2vWrFl68sknK3VWAMCNrax7X1o1vFWdwoP1/Mo9MjLqGB6s1MyTmrryW93TrI5aNbz1+g9qoRs+aGrVungKr6xLUGlpaaWWDR48WMePHy+1fMeOHUpLS9O9996rJk2aKCIiovKHBQDclOYOa68Ji9M0PunfP5fuaVZHsweXvu0CZbvhg6ZDh4t3h7/33nsaPny484zNli1b9MEHH5Tavry3fCckJCgtLU2PPfaYBg8eXHUDAwBuOsH+Plo4qpMyjmYr60SuGteuoWahAe4eyyo3fNB07txZUVFR+vzzzxUVFaWuXbvq+++/16pVq9SvXz+tWLHC3SMCACBJahYaQMhcpWr/Lqdr5XA4tGrVKg0fPlwHDhzQn//8Z/3www9atWoV72ICAOAG4TDGGHcPAQAAcC1u+DM0AADgxkfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKx3XYLmzTff1G9+8xu1atVKXl5ecjgc2rBhw/V4aJcZY7RmzRo9/vjjatWqlYKCglSjRg21bt1aM2bMUF5enrtHBAAAv+AwxpgqfxCHQ5IUFhYmY4x++uknrV+/XrGxsVX90C7Ly8vTLbfcIl9fX8XGxuruu+9WXl6eUlJSlJGRocjISG3cuFG33HKLu0cFAAD/57qcofn444/1z3/+U0eOHFH//v2vx0NeNU9PT02fPl0//fST1q5dq5kzZ+rNN9/Ut99+q379+mnbtm3605/+5O4xr5uMo9n6296jyjia7e5RAAAol8tBU1hYqJkzZ6pFixby8/NTkyZNNG3aNBUWFmrv3r1yOByaN29eiX369OmjevXqVdrQrjp06JCGDx+uOnXqKCAgQNHR0Vq3bp0k6emnn5afn5+ysy/+wPb29tazzz6rW2+9tcQxvL299cc//lGStHHjxus6//VyIiff+XXgWLYGvbVFPWdt0m/f+1o9Z23SoLe26MCxbOc2AABUF16u7jBo0CAtW7ZM0dHRio+P18aNG5WQkKCCggJ5eXnJw8ND8fHxVTDq1Tl48KCio6N16tQpDR06VEFBQVq4cKH69u2r1NRULV26VD169FBAQMAVj+Xt7S1J8vJy+a/NrXLPX6jQdu1f+sz5Zw+HFODnrT8PaafI8FralnlKz674Rr1mbVLR/12k3PtiXIWOW8PHrr8vAIB9XPpJs3z5ci1btkxdu3bV+vXr5eHhoaKiIkVFRWnu3LkKCQlRTEyM6tatW1XzumzixIk6duyY3n77bY0aNUqSNHLkSLVt21Zjx45VVlaWnn/++Qoda/78+ZKkXr16Vdm8VeHOqSku71NkpBkD7lafVmGSpD6twmRkND4pzeXjZr3Sx+XHBwDAFS5dclqyZImki5Hg4XFxVw8PD40ePVonT57Ud999pwceeKDyp7xKubm5+uSTTxQSEqLhw4c7l7dq1UqdO3fW5s2b5enpWaH7etauXau33npLLVu2dIbRjS4yvFaJ7zuGB7tpEgAALs+lMzT79u2TJMXExJRYHhUV5fzzgAEDKmGs0pKTk7Vz584Sy2JjYy/7Tqn9+/ersLBQkZGR8vHxKbEuKipKmzdvVrdu3VS7du3LPvbXX3+tQYMGKSgoSEuWLJGvr+/VPg23qOiloV+ecdmWecp5hkaSUjNPXtVxAQCoai4FTU5Ojjw9PRUSElJiecuWLeXn56eIiAg1atSoUgcslpycrHfffbfU8ssFTU5OjqSLbxf/pbZt20qSBg4ceNnHTUtLU69eveRwOJSSkqK77rrLhamrh4rew7L9uR7OP4/7YIeeX7lHRkYdw4OVmnlSU1d+q07hwZoztJ1LxwUAoKq59BPJ399fhYWFysvLk5+fn3N5Zmam8vLy1KBBg0ofsFhiYqISExNd2sff31+S9PPPP5daV3y26XIz79ixQz179lRhYaE+/fRTRUZGuvT4tqld899nnuYOa68Ji9NK3DNzT7M6mj24rYL9fcraHQAAt3EpaFq2bKndu3drx44dio6Odi5fuHChJGnXrl2VO901at68uTw9PbV9+/YSy4uKipSUlCTp4sxlvStrx44d6tGjhy5cuKCUlBR16tTpeoxcbQT7+2jhqE7KOJqtrBO5aly7hpqFXvmdYAAAuINLNwUX/+BPSEjQ+fPnJUknT57UW2+9JUnKysrSd999V7kTXgN/f3/16NFDGRkZJc7uLFq0SJmZmZKklJTS79QpjpmCggKtWbOmxD1CN5tmoQHqeWcoMQMAqNZc+tUHRUVF6tq1q7788ku1bt1a3bt315o1a5Senq4ZM2ZoypQpioiI0JAhQ/TMM88493vllVeUnp4uSdqyZYv279+vuLg454ftjR49Wl26dKnkp3ZRWlqaYmJilJ+frwceeEAhISGaP3++QkJC1K9fP82ZM0fDhg3T0KFD1bt3b508eVJNmzbVqVOn1Lt37zLPzNx666168sknq2ReAABwFYyLTp8+bSZMmGAaNmxovL29TWhoqJk9e7Yxxph58+aZ0NBQExQUVGKfbt26GUnlfi1YsMDVMVyydetWExcXZwIDA42Pj49p06aN2bNnjzl79qwZPHiw8fX1NRMmTDDGGJOZmXnZWSWZ22+/vUrnBQAArrkuv5wSAACgKl2XX04JAABQlQgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPWuS9C8+eab+s1vfqNWrVrJy8tLDodDGzZsuB4PfVWWL1+uBx98UM2aNVNgYKBq1qypu+66S08++aQOHz7s7vEAAMAvOIwxpsofxOGQJIWFhckYo59++knr169XbGxsVT/0VRk+fLi++uorRUZGKiwsTJK0c+dOrV+/XkFBQfriiy901113uXlKAABQ7LoEzerVq9W+fXvVq1dPY8eO1VtvvVWtgyYvL09+fn6llr/zzjsaPXq0HnzwQS1ZssQNk1WOjKPZyjqRq8a1a6hZaIC7xwEA4Jq5fMmpsLBQM2fOVIsWLeTn56cmTZpo2rRpKiws1N69e+VwODRv3rwS+/Tp00f16tWrtKFddejQIQ0fPlx16tRRQECAoqOjtW7dOknS008/LT8/P2VnZzu3LytmJOnXv/61JOnAgQNVP3QlOpGTrxM5+TpwLFuD3tqinrM26bfvfa2eszZp0FtbdOBYtnMbAABs5OXqDoMGDdKyZcsUHR2t+Ph4bdy4UQkJCSooKJCXl5c8PDwUHx9fBaNenYMHDyo6OlqnTp3S0KFDFRQUpIULF6pv375KTU3V0qVL1aNHDwUEXPlMxerVqyVJERERVT12heSev1Ch7dq/9JkkycMhBfh5689D2ikyvJa2ZZ7Ssyu+Ua9Zm1T0f+fp9r4YV+HHr+Hj8v98AACoEi79RFq+fLmWLVumrl27av369fLw8FBRUZGioqI0d+5chYSEKCYmRnXr1q2qeV02ceJEHTt2TG+//bZGjRolSRo5cqTatm2rsWPHKisrS88//3yZ+yYnJ2vnzp3Kzc3Vt99+q5SUFIWHh+vFF1+8nk+hXHdOTXFp+yIjzRhwt/q0unhfUJ9WYTIyGp+UdlXHzHqlj0uPDwBAVXEpaIrvG5k4caI8PC5erfLw8NDo0aP12GOP6eTJk3r88ccrf8qrlJubq08++UQhISEaPny4c3mrVq3UuXNnbd68WZ6enurfv3+Z+ycnJ+vdd991ft+hQwctXrxY4eHhVT57VYkMr1Xi+47hwW6aBACAyuNS0Ozbt0+SFBMTU2J5VFSU888DBgyohLFKKz5bcqnY2NjL3li8f/9+FRYWKjIyUj4+PiXWRUVFafPmzerWrZtq165d5v6JiYlKTEzUmTNnlJaWpilTpqh9+/Zavny5unfvfq1P6ZpV9PLQpWddtmWecp6hkaTUzJNXdUwAAKoTl4ImJydHnp6eCgkJKbG8ZcuW8vPzU0REhBo1alSpAxb75dmSYpcLmpycHElyvvX6Um3btpUkDRw48IqPHRQUpNjYWK1Zs0YtWrTQI488oszMTHl7e1dw+qpR0XtYtj/XQ5I07oMden7lHhkZdQwPVmrmSU1d+a06hQdrztB2Lh0TAIDqxKWfXv7+/iosLCz1tubMzEzl5eWpQYMGlT5gseKzJa7w9/eXJP3888+l1hWfbXJl5sDAQHXu3FnJyck6cOCAWrZs6dI87lK7pq8kae6w9pqwOK3EPTP3NKuj2YPbKtjfp7zdAQCo9lwKmpYtW2r37t3asWOHoqOjncsXLlwoSdq1a1flTneNmjdvLk9PT23fvr3E8qKiIiUlJUm6OLMr78o6cuSIJMnLy74zGcH+Plo4qhOfQwMAuOG49Dk0xT/4ExISdP78eUnSyZMn9dZbb0mSsrKy9N1331XuhNfA399fPXr0UEZGRomzO4sWLVJmZqYkKSWl5Lt68vPz9dVXX5V5vAULFig1NVVNmzZVs2bNqmzuqtYsNEA97wwlZgAANwyXPim4qKhIXbt21ZdffqnWrVure/fuWrNmjdLT0zVjxgxNmTJFERERGjJkiJ555hnnfq+88orS09MlSVu2bNH+/fsVFxfn/LC90aNHq0uXLpX81C5KS0tTTEyM8vPz9cADDygkJETz589XSEiI+vXrpzlz5mjYsGEaOnSoevfurdOnT6tWrVqKiIhQmzZt1KBBA505c0apqanasWOHatasqTVr1lTZvAAA4CoYF50+fdpMmDDBNGzY0Hh7e5vQ0FAze/ZsY4wx8+bNM6GhoSYoKKjEPt26dTOSyv1asGCBq2O4ZOvWrSYuLs4EBgYaHx8f06ZNG7Nnzx5z9uxZM3jwYOPr62smTJhgjDHm/PnzZtq0aSY2NtaEhYUZb29vU6NGDXPnnXeaJ5980nz//fdVOisAAHDddfldTgAAAFXJ5d/lBAAAUN0QNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsd12C5s0339RvfvMbtWrVSl5eXnI4HNqwYcP1eOhKcfr0aTVo0EAOh0O9e/d29zgAAOAXvK7HgzzxxBOSpLCwMIWEhOinn366Hg9baZ544gmdOXPG3WMAAIByXJczNB9//LH++c9/6siRI+rfv//1eMhK89FHH2nhwoV6+eWX3T0KAAAoh8tBU1hYqJkzZ6pFixby8/NTkyZNNG3aNBUWFmrv3r1yOByaN29eiX369OmjevXqVdrQrjp06JCGDx+uOnXqKCAgQNHR0Vq3bp0k6emnn5afn5+ys7NL7Xfy5Ek99thjGjJkiPr163e9x76pZBzN1t/2HlXG0dKvAwAAV+LyJadBgwZp2bJlio6OVnx8vDZu3KiEhAQVFBTIy8tLHh4eio+Pr4JRr87BgwcVHR2tU6dOaejQoQoKCtLChQvVt29fpaamaunSperRo4cCAgJK7Tt+/HgVFhbqjTfeKDN4cPVO5ORLkk7lnteUFXu0NfOkc12n8GBNHxChWjV8JEm1a/q6ZUYAgD1cCprly5dr2bJl6tq1q9avXy8PDw8VFRUpKipKc+fOVUhIiGJiYlS3bt2qmtdlEydO1LFjx/T2229r1KhRkqSRI0eqbdu2Gjt2rLKysvT888+X2m/FihVatGiRPvzwQ9WuXZugqaDc8xcqtF37lz6TJHk4pAA/b/15SDtFhtfStsxTenbFN+o1a5OKzMVt974YV6Fj1vC5LreEAQCqIZd+AixZskTSxUjw8Lh4tcrDw0OjR4/WY489ppMnT+rxxx+v/CmvUm5urj755BOFhIRo+PDhzuWtWrVS586dtXnzZnl6epa6r+f48eMaO3as4uPj9dBDD13vsa1259QUl7YvMtKMAXerT6swSVKfVmEyMhqflObyMbNe6ePSYwMAbhwuBc2+ffskSTExMSWWR0VFOf88YMCAShirtOTkZO3cubPEstjYWMXGxpa7z/79+1VYWKjIyEj5+PiUWBcVFaXNmzerW7duql27dol148aNU0FBgebOnVtZ4+MyIsNrlfi+Y3iwmyYBANjKpaDJycmRp6enQkJCSixv2bKl/Pz8FBERoUaNGlXqgMWSk5P17rvvllp+uaDJycmRdPHt4r/Utm1bSdLAgQNLLF+5cqWWLFmixMREt97IbKuKXh669KzLtsxTzjM0kpR6yf00rhwTAHDzcilo/P39VVhYqLy8PPn5+TmXZ2ZmKi8vTw0aNKj0AYslJiYqMTHRpX38/f0lST///HOpdcVnm345c1raxUsdI0aM0IgRI0rtl5KSIofDodatW5c6Y4SK38ey/bkekqRxH+zQ8yv3yMioY3iwUjNPaurKb9UpPFhzhrZz6ZgAgJuXSz8pWrZsqd27d2vHjh2Kjo52Ll+4cKEkadeuXZU73TVq3ry5PD09tX379hLLi4qKlJSUJOnizJe+K6tdu3bOm4cvlZOTow8//FANGzZUXFxclZ2JulkUv3Np7rD2mrA4rcQ9M/c0q6PZg9sq2N+nvN0BACjJuGDRokVGkunZs6fJz883xhhz4sQJExoaaiQZSSY9Pf2yxxgzZoyRZNavX+/KQ1+1uLg4I8ksWLDAuez99993zhsVFVWh42RmZhpJJi4uroomvbnt/+ms+fTbn8z+n866exQAgIUcxhhT0fgpKipS165d9eWXX6p169bq3r271qxZo/T0dM2YMUNTpkxRRESEhgwZomeeeca53yuvvKL09HRJ0pYtW7R//37FxcU571EZPXq0unTpUkmJVlJaWppiYmKUn5+vBx54QCEhIZo/f75CQkLUr18/zZkzR8OGDdPQoUMv+3uasrKyFB4erri4OK1du7ZKZgUAAFfJ1QI6ffq0mTBhgmnYsKHx9vY2oaGhZvbs2cYYY+bNm2dCQ0NNUFBQiX26devmPCNS1telZ0+qwtatW01cXJwJDAw0Pj4+pk2bNmbPnj3m7NmzZvDgwcbX19dMmDDhssfgDA0AANWXS2doAAAAqqPr8sspAQAAqhJBAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALBelQdNRkaGZsyYoa5du6p+/fry8fHRbbfdpkceeUTp6elV/fBXZfny5XrwwQfVrFkzBQYGqmbNmrrrrrv05JNP6vDhw+4eDwAA/ILDGGOq8gEGDx6sDz/8UBEREerSpYsCAwP1zTffaM2aNbrllluUkpKie+65pypHcNnw4cP11VdfKTIyUmFhYZKknTt3av369QoKCtIXX3yhu+66y81TAgCAYlUeNImJiWrbtq1at25dYvnixYv18MMP684779S3335blSO4LC8vT35+fqWWv/POOxo9erQefPBBLVmyxA2TAQCAsrh8yamwsFAzZ85UixYt5OfnpyZNmmjatGkqLCzU3r175XA4NG/ePOf2I0aMKBUz0sUzN82bN9fevXt1/Pjxa3sWV3Do0CENHz5cderUUUBAgKKjo7Vu3TpJ0tNPPy0/Pz9lZ2c7ty8rZiTp17/+tSTpwIEDVTrvjSTjaLb+tveoMo5mX3ljAACukperOwwaNEjLli1TdHS04uPjtXHjRiUkJKigoEBeXl7y8PBQfHx8hY7l7e19cQgvl8eosIMHDyo6OlqnTp3S0KFDFRQUpIULF6pv375KTU3V0qVL1aNHDwUEBFzxWKtXr5YkRUREVNm8tjuRky9JOpV7XlNW7NHWzJPOdZ3CgzV9QIRq1fBR7Zq+7hoRAHADcumS0/Lly/XAAw+oa9euWr9+vTw8PFRUVKSoqCgdOHBAISEhqlu3rjZt2nTFY6WmpqpTp06KjIxUamrqNT2Jy7n//vv10Ucf6e2339aoUaMkSbt371bbtm3VuXNnbd68We+8845GjhxZat/k5GTt3LlTubm5+vbbb5WSkqJGjRpp3bp1Cg8Pr7KZq6Pc8xcqtN2dU1MkSR4OKcDPWzMG3K3I8FralnlKz674Rtl5BSoy0t4X4654rBo+VRe6AIAbi0tB8/DDD2vx4sVasWJFibMw//u//6vHHntMkvQ///M/mjBhwmWPc+bMGXXu3Fn79+/XunXrFBsbe1XDX0lubq4CAwMVHBysH3/8UT4+Ps51MTEx2rx5szw9PXX06FHVrl271P4jRozQu+++6/y+Q4cOWrx4se64444qmbc6a/zMapf3+fOQdurTKsz5/ce7j2h8UlqF9896pY/LjwkAuDm59E/gffv2SboYA5eKiopy/nnAgAGXPUZeXp4GDhyo9PR0TZ8+vcIxU3y25FKxsbGX3X///v0qLCxUZGRkiZgpnnnz5s3q1q1bmTEjXbyhOTExUWfOnFFaWpqmTJmi9u3ba/ny5erevXuF5r6ZRYbXKvF9x/BgN00CALjRuRQ0OTk58vT0VEhISInlLVu2lJ+fnyIiItSoUaNy98/Pz9eAAQP0+eef649//KOeffbZCj92cnJyibMlxS4XNDk5OZLkfOv1pdq2bStJGjhw4BUfOygoSLGxsVqzZo1atGihRx55RJmZmc57gG4GFblEJP37kpMkbcs8VeIMTeol99NU9HgAAFSES0Hj7++vwsLCUm9rzszMVF5enho0aFDuvnl5eYqPj1dKSor+8Ic/aMaMGS4NWny2xNV5Jennn38uta74bNPlZv6lwMBAde7cWcnJyTpw4IBatmzp0jw2q+j9LNuf6yFJGvfBDj2/co+MjDqGBys186SmrvxWncKDNWdoO+6PAQBUKpfetl38A3zHjh0lli9cuFCStGvXrjL3uzRmnnrqKb366qtXM6vLmjdvLk9PT23fvr3E8qKiIiUlJUkqf+byHDlyRFLVvjPLZrVr+qp2TV/NHdZed9UP1PikNHWcvk7jk9J0V/1AzR3Wnnc4AQAqn3HBokWLjCTTs2dPk5+fb4wx5sSJEyY0NNRIMpJMenp6iX3OnTtnevXqZSSZSZMmufJwlSIuLs5IMgsWLHAue//9953zRkVFldg+Ly/PbNmypcxjzZ8/30gyTZs2rcqRbyj7fzprPv32J7P/p7PuHgUAcANz6V1ORUVF6tq1q7788ku1bt1a3bt315o1a5Senq4ZM2ZoypQpioiI0JAhQ/TMM89I+vc7herVq6cxY8aUedwRI0aocePG1xRm5UlLS1NMTIzy8/P1wAMPKCQkRPPnz1dISIj69eunOXPmaNiwYRo6dKh69+6t06dPq1atWoqIiFCbNm3UoEEDnTlzRqmpqdqxY4dq1qypNWvWqEuXLlUyLwAAuAquFtDp06fNhAkTTMOGDY23t7cJDQ01s2fPNsYYM2/ePBMaGmqCgoKc23fr1s15NqS8r/Xr11dOnpVj69atJi4uzgQGBhofHx/Tpk0bs2fPHnP27FkzePBg4+vrayZMmGCMMeb8+fNm2rRpJjY21oSFhRlvb29To0YNc+edd5onn3zSfP/991U6KwAAcF2V/y4nAACAquby73ICAACobggaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPWqPGgyMjI0Y8YMde3aVfXr15ePj49uu+02PfLII0pPT6/qh68Up0+fVoMGDeRwONS7d293jwMAAH7Bq6of4Pnnn9eHH36oiIgI9e/fX4GBgfrmm2+0cOFCLV26VCkpKbrnnnuqeoxr8sQTT+jMmTPuHgMAAJSjyoOmd+/e+uMf/6jWrVuXWL548WI9/PDDGjt2rL799tuqHuOqffTRR1q4cKHeeOMNPfHEE+4eBwAAlMHlS06FhYWaOXOmWrRoIT8/PzVp0kTTpk1TYWGh9u7dK4fDoXnz5jm3HzFiRKmYkaTBgwerefPm2rt3r44fP35tz+IKDh06pOHDh6tOnToKCAhQdHS01q1bJ0l6+umn5efnp+zs7FL7nTx5Uo899piGDBmifv36VemMuH4yjmbrb3uPKuNo6dccAGAnl8/QDBo0SMuWLVN0dLTi4+O1ceNGJSQkqKCgQF5eXvLw8FB8fHyFjuXt7X1xCK+qO1F08OBBRUdH69SpUxo6dKiCgoK0cOFC9e3bV6mpqVq6dKl69OihgICAUvuOHz9ehYWFeuONN8oMHtjhRE6+JOlU7nlNWbFHWzNPOtd1Cg/W9AERqlXDR7Vr+rprRADANXKpJJYvX65ly5apa9euWr9+vTw8PFRUVKSoqCjNnTtXISEhiomJUd26da94rNTUVH377beKjIzUrbfeerXzX9HEiRN17Ngxvf322xo1apQkaeTIkWrbtq3Gjh2rrKwsPf/886X2W7FihRYtWqQPP/xQtWvXJmiqqdzzF664TfuXPpMkeTikAD9v/XlIO0WG19K2zFN6dsU36jVrk4qMtPfFuMsep4ZPlV+hBQBcJZf+C71kyRJJFyPBw+Pi1SoPDw+NHj1ajz32mE6ePKnHH3/8isc5c+aMHn30UXl4eOi11167irErJjc3V5988olCQkI0fPhw5/JWrVqpc+fO2rx5szw9PdW/f/8S+x0/flxjx45VfHy8HnrooSqbD9fuzqkpFd62yEgzBtytPq3CJEl9WoXJyGh8UlqFjpX1Sp+rHxQAUKVcCpp9+/ZJkmJiYkosj4qKcv55wIABlz1GXl6eBg4cqPT0dE2fPl2xsbEVeuzk5GTt3LmzxLLY2NjL7r9//34VFhYqMjJSPj4+pWbevHmzunXrptq1a5dYN27cOBUUFGju3LkVmg32iAyvVeL7juHBbpoEAFCZXAqanJwceXp6KiQkpMTyli1bys/PTxEREWrUqFG5++fn52vAgAH6/PPP9cc//lHPPvtshR87OTlZ7777bqnllwuanJwcSVJYWFipdW3btpUkDRw4sMTylStXasmSJUpMTFS9evUqPB/c40qXiaSSZ162ZZ5ynqGRpNRL7qepyLEAANWTS0Hj7++vwsJC5eXlyc/Pz7k8MzNTeXl5atCgQbn75uXlKT4+XikpKfrDH/6gGTNmuDRoYmKiEhMTXdrH399fkvTzzz+XWld8tumXM6elXbz8MGLECI0YMaLUfikpKXI4HGrdunWpM0a4/ipyX8v253pIksZ9sEPPr9wjI6OO4cFKzTypqSu/VafwYM0Z2o57ZADAYi79F7xly5bavXu3duzYoejoaOfyhQsXSpJ27dpV5n6XxsxTTz2lV1999RpGrrjmzZvL09NT27dvL7G8qKhISUlJki7OfOm7stq1a+e8efhSOTk5+vDDD9WwYUPFxcVd9kwUqpfidy/NHdZeExanOe+ZkaR7mtXR7MFtFezvU97uAAAbGBcsWrTISDI9e/Y0+fn5xhhjTpw4YUJDQ40kI8mkp6eX2OfcuXOmV69eRpKZNGmSKw9XKeLi4owks2DBAuey999/3zlvVFRUhY6TmZlpJJm4uLgqmhTXy/6fzppPv/3J7P/prLtHAQBUEocxxlQ0foqKitS1a1d9+eWXat26tbp37641a9YoPT1dM2bM0JQpUxQREaEhQ4bomWeekXTx0s27776revXqacyYMWUed8SIEWrcuPE1hVl50tLSFBMTo/z8fD3wwAMKCQnR/PnzFRISon79+mnOnDkaNmyYhg4detnf05SVlaXw8HDFxcVp7dq1VTIrAAC4Sq4W0OnTp82ECRNMw4YNjbe3twkNDTWzZ882xhgzb948ExoaaoKCgpzbd+vWzXk2pLyv9evXV06elWPr1q0mLi7OBAYGGh8fH9OmTRuzZ88ec/bsWTN48GDj6+trJkyYcNljcIYGAIDqy6UzNAAAANWRy7/LCQAAoLohaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWq/Kg+f777zV27Fi1b99eISEh8vX11e23364+ffpo3bp1Vf3wleL06dNq0KCBHA6Hevfu7e5xAADAL1R50GRkZOjDDz9UaGioHnroIU2aNEndunXTF198oR49emjGjBlVPcI1e+KJJ3TmzBl3jwEAAMrhMMaYqnyA8+fPy8vLSx4eJdvpyJEjateunU6ePKljx47p1ltvrcoxrtpHH32k+++/X2+88YaeeOIJxcXFae3ate4eC9dJxtFsZZ3IVePaNdQsNMDd4wAAyuHyGZrCwkLNnDlTLVq0kJ+fn5o0aaJp06apsLBQe/fulcPh0Lx585zb+/j4lIoZSapfv76io6NVUFCg77///tqexRUcOnRIw4cPV506dRQQEKDo6Gjn5a6nn35afn5+ys7OLrXfyZMn9dhjj2nIkCHq169flc6I6uFETr5O5OTrwLFsDXpri3rO2qTfvve1es7apEFvbdGBY6X/dwIAcD8vV3cYNGiQli1bpujoaMXHx2vjxo1KSEhQQUGB80xMfHz8FY9z4sQJbd26VTVq1FCTJk2uZvYKOXjwoKKjo3Xq1CkNHTpUQUFBWrhwofr27avU1FQtXbpUPXr0UEBA6X99jx8/XoWFhXrjjTfKDB7YJff8hStu0/6lzyRJHg4pwM9bfx7STpHhtbQt85SeXfGNes3apD3T4q54nBo+Lv9fCwBwDVz6r+7y5cu1bNkyde3aVevXr5eHh4eKiooUFRWluXPnKiQkRDExMapbt26pfbOyspSYmKjCwkIdOXJEq1at0unTp/WXv/ylzJioLBMnTtSxY8f09ttva9SoUZKkkSNHqm3btho7dqyysrL0/PPPl9pvxYoVWrRokT788EPVrl2boLkB3Dk1pcLbFhlpxoC71adVmCSpT6swGRmNT0qr0HGyXulz1XMCAFznUtAsWbJE0sVIKL6M5OHhodGjR+uxxx7TyZMn9fjjj5e5b1ZWlqZNm+b8vmbNmlqwYIGGDRt2tbNfUW5urj755BOFhIRo+PDhzuWtWrVS586dtXnzZnl6eqp///4l9jt+/LjGjh2r+Ph4PfTQQ1U2H6q3yPBaJb7vGB7spkkAAFfiUtDs27dPkhQTE1NieVRUlPPPAwYMKHPf2NhYGWNUUFCgrKwszZs3T4888ohSU1P1xhtvXPGxk5OTtXPnzlLHjI2NLXef/fv3q7CwUJGRkfLx8Sk18+bNm9WtWzfVrl27xLpx48apoKBAc+fOveJcsMfeF698qejSsy/bMk85z9BIUmrmSUnSqvExalq3ZuUPCAC4ai4FTU5Ojjw9PRUSElJiecuWLeXn56eIiAg1atTossfw9vZWs2bNNHPmTOXm5urNN9/Ufffdp/vuu++y+yUnJ+vdd98ttfxyQZOTkyNJCgsLK7Wubdu2kqSBAweWWL5y5UotWbJEiYmJqlev3mVngl0qcl/L9ud6SJLGfbBDz6/cIyOjjuHBSs08qakrv1Wn8GC1anhrFU8KAHCVS+9y8vf3V2FhofLy8kosz8zMVF5enho0aODSg/fq1UuStGHDhitum5iYKGNMia+EhIQrzitJP//8c6l1xWebfjlzWlqaJGnEiBFyOBzOr/DwcElSSkqKHA6H2rRpc8WZYZ/aNX1Vu6av5g5rr7vqB2p8Upo6Tl+n8Ulpuqt+oOYOa+/uEQEAZXDpDE3Lli21e/du7dixQ9HR0c7lCxculCTt2rXLpQc/cuTIxSG8quYdIc2bN5enp6e2b99eYnlRUZGSkpIkXZz50ndltWvXznnz8KVycnL04YcfqmHDhoqLi7vimSjYLdjfRwtHdeJzaADAFsYFixYtMpJMz549TX5+vjHGmBMnTpjQ0FAjyUgy6enpJfbZunWrOXfuXKljZWVlmdtuu81IMn//+99dGcMlcXFxRpJZsGCBc9n777/vnDcqKqpCx8nMzDSSTFxcXBVNCgAArpZLnxRcVFSkrl276ssvv1Tr1q3VvXt3rVmzRunp6ZoxY4amTJmiiIgIDRkyRM8884wkKT4+Xn//+9/VrVs3NWrUSF5eXjp48KA++eQTnT9/XhMnTtR///d/V3qoFUtLS1NMTIzy8/P1wAMPKCQkRPPnz1dISIj69eunOXPmaNiwYRo6dOhlf09TVlaWwsPD+aRgAACqI1cL6PTp02bChAmmYcOGxtvb24SGhprZs2cbY4yZN2+eCQ0NNUFBQc7tP/roIzN48GBzxx13GH9/f+Pt7W0aNGhgBgwYYD755JNK6rLL27p1q4mLizOBgYHGx8fHtGnTxuzZs8ecPXvWDB482Pj6+poJEyZc9hicoQEAoPqq8t/lBAAAUNWq/LdtAwAAVDWCBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWI+gAQAA1iNoAACA9QgaAABgPYIGAABYj6ABAADWI2gAAID1CBoAAGA9ggYAAFiPoAEAANYjaAAAgPUIGgAAYD2CBgAAWM/L3QPcDIwxys7OdvcYAABYKyAgQA6Ho9z1BM11kJ2draCgIHePAQCAtc6cOaPAwMBy1zuMMeY6znNT4gxNxZw9e1a33Xabfvjhh8v+jxbuw2tkB14nO/A6uYYzNNWAw+Hgf6wuCAwM5O+rmuM1sgOvkx14nSoHNwUDAADrETQAAMB6BA2qDV9fX73wwgvy9fV19ygoB6+RHXid7MDrVLm4KRgAAFiPMzQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNHC7bdu26T/+4z9Uq1Yt+fv7q2PHjkpKSnL3WDeVw4cP63/+53/Uq1cvNWrUSD4+PqpXr54eeOABbd26tcx9zp49q0mTJun222+Xr6+vbr/9dk2aNElnz569ztPf3F577TU5HA45HA599dVXZW7Da+U+K1asUM+ePVW7dm3dcsstCg8P18MPP6wffvihxHa8RteOdznBrTZs2KC4uDj5+Pho8ODBCgoK0vLly5WZmanp06fr2WefdfeIN4VnnnlGr776qu644w5169ZNdevWVUZGhpKTk2WM0aJFi/TQQw85t//555/VpUsX7dy5Uz179lS7du20a9curV27Vm3atNEXX3whf39/Nz6jm8O+ffvUtm1beXl56eeff9aWLVvUuXPnEtvwWrmHMUZjx47VvHnzdMcddyguLk4BAQE6cuSINm7cqA8++EBdunSRxGtUaQzgJgUFBeaOO+4wvr6+ZseOHc7lZ8+eNXfddZfx8vIy+/fvd+OEN49ly5aZTZs2lVq+adMm4+3tbYKDg01eXp5z+dSpU40k84c//KHE9sXLp06dWuUz3+wuXLhgIiMjTceOHc2wYcOMJLNly5ZS2/Faucfs2bONJPO73/3OXLhwodT6goIC5595jSoHQQO3SUlJMZLMb37zm1LrFi9ebCSZP/7xj26YDJfq1auXkWS2bdtmjDGmqKjI1K9f39SsWdPk5OSU2PbcuXOmVq1apkGDBqaoqMgd4940pk+fbnx8fMyePXvMo48+WmbQ8Fq5R25urgkODjZNmjQpES5l4TWqPNxDA7fZsGGDJKlXr16l1hUv27hx4/UcCWXw9vaWJHl5XfxdthkZGTpy5IhiYmJKnQb38/NT165ddfjwYR04cOC6z3qz2LNnj6ZNm6bnnntOd911V7nb8Vq5x9/+9jedPHlS8fHxKiws1PLly/XKK6/oL3/5S6m/a16jysNv24bbZGRkSJKaNWtWal2tWrVUp04d5zZwj0OHDumzzz5TvXr1dPfdd0u6/Ot26fKMjIxyt8HVu3DhgkaMGKGWLVvqmWeeuey2vFbu8fXXX0u6+I+A1q1b67vvvnOu8/Dw0MSJE/Vf//VfkniNKhNnaOA2Z86ckSQFBQWVuT4wMNC5Da6/goICDR8+XPn5+Xrttdfk6ekpqWKv26XboXLNmDFDu3bt0vz5851nz8rDa+Uex44dkyS9/vrrCgwMVGpqqrKzs7Vp0yY1b95cr7/+uubOnSuJ16gyETQASikqKtLIkSO1adMm/fa3v9Xw4cPdPRIk7dq1Sy+99JKeeuoptWvXzt3joBxFRUWSJB8fHyUnJysyMlI1a9bUPffco6VLl8rDw0Ovv/66m6e88RA0cJvif5GU9y+Ps2fPlvuvFlQdY4x++9vf6v3339ewYcP0l7/8pcT6irxul26HyvPoo4/qjjvuUEJCQoW257Vyj+K/zw4dOqh+/fol1t11111q0qSJDh48qNOnT/MaVSKCBm5z6bXhXzp16pSOHz/ONePrrKioSKNGjdL8+fP18MMPKzExUR4eJf8zcbnX7dLlvHaVb9euXUpPT5efn5/zw/QcDofeffddSVJUVJQcDoeSk5Ml8Vq5S4sWLSRJt956a5nri5efO3eO16gScVMw3KZbt256+eWX9emnn2rw4MEl1n366afObXB9FBUVafTo0VqwYIEGDRqkhQsXOu+buVSzZs1Uv359ffnll/r5559LvDMjLy9PmzZtUv369dW0adPrOf5NYdSoUWUu37RpkzIyMnT//fcrJCREjRs3lsRr5S733nuvpIsffPhLBQUFOnDggPz9/RUSEqJ69erxGlUWd79vHDevgoIC06RJE+Pr62vS0tKcyy/9YL3vvvvOfQPeRAoLC82IESOMJPPrX//6ip+dwQeBVS/lfQ6NMbxW7lL8+U3/+7//W2L5iy++aCSZYcOGOZfxGlUOfvUB3Gr9+vWKi4uTr6+vHn74YQUGBjp/9cFLL72kKVOmuHvEm0JCQoKmTZummjVrasKECc7PnLlUfHy82rRpI6n0R7W3b99eu3bt0po1a/iodjcYMWKE3n333Qr96gNeq+vj4MGDio6O1rFjx9SnTx/9v//3/5SWlqbPP/9ct99+u7766ivVq1dPEq9RpXF3UQFbt241vXv3NkFBQeaWW24xHTp0MO+//767x7qpFP8L/3JfCxYsKLHP6dOnzcSJE81tt91mvL29zW233WYmTpxoTp8+7Z4ncRO73BkaY3it3OXQoUNmxIgRpl69es6/99/97nfm6NGjpbblNbp2nKEBAADW411OAADAegQNAACwHkEDAACsR9AAAADrETQAAMB6BA0AALAeQQMAAKxH0AAAAOsRNAAAwHoEDQAAsB5BAwAArEfQAAAA6xE0AADAev8/2AQgSHChkGgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "az.plot_forest(\n", " trace, var_names=[\"~τ\", \"~μ0\", \"~α1\", \"~α2\", \"~α3\", \"~α4\"], combined=True\n", ")" ] }, { "cell_type": "markdown", "id": "ec7443cc", "metadata": {}, "source": [ "## A more concise method\n", "Not necessarily pretty, but easier to extend to more treatments. I'm interested in seeing other peoples' methods here; I feel like this could still be a lot cleaner." ] }, { "cell_type": "code", "execution_count": 7, "id": "3b43b405", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# get possible combinations\n", "combos = list(combinations(range(4), 2))\n", "combos" ] }, { "cell_type": "code", "execution_count": 8, "id": "9801d470", "metadata": { "tags": [ "hide-output" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [μ0, τ, α2, α3, α4]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [24000/24000 00:03<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 4 seconds.\n" ] } ], "source": [ "with pm.Model() as m:\n", " mu0 = pm.Normal(\"μ0\", mu=0, tau=0.0001)\n", " tau = pm.Gamma(\"τ\", 0.001, 0.001)\n", "\n", " alphas = [pm.Normal(f\"α{i}\", mu=0, tau=0.0001) for i in range(2, 5)]\n", "\n", " # sum-to-zero constraint\n", " alphas.insert(0, pm.Deterministic(\"α1\", -(alphas[0] + alphas[1] + alphas[2])))\n", "\n", " mus = [\n", " pm.Deterministic(f\"mu{i + 1}\", mu0 + alpha) for i, alpha in enumerate(alphas)\n", " ]\n", "\n", " likelihoods = [\n", " pm.Normal(f\"lik{i + 1}\", mu=mus[i], tau=tau, observed=data[i + 1])\n", " for i, mu in enumerate(mus)\n", " ]\n", "\n", " [pm.Deterministic(f\"α{i + 1} - α{j + 1}\", alphas[i] - alphas[j]) for i, j in combos]\n", "\n", " trace = pm.sample(5000)" ] }, { "cell_type": "code", "execution_count": 9, "id": "f8e34c05", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_3%hdi_97%
α22.0040.8910.3203.656
α34.0090.8882.3025.631
α4-3.0030.805-4.530-1.489
α1-3.0101.035-4.925-0.993
α1 - α2-5.0141.621-8.095-2.047
α1 - α3-7.0191.617-9.971-3.907
α1 - α4-0.0071.532-2.9142.883
α2 - α3-2.0051.438-4.7750.659
α2 - α45.0071.3342.4317.451
α3 - α47.0111.3314.4899.524
\n", "
" ], "text/plain": [ " mean sd hdi_3% hdi_97%\n", "α2 2.004 0.891 0.320 3.656\n", "α3 4.009 0.888 2.302 5.631\n", "α4 -3.003 0.805 -4.530 -1.489\n", "α1 -3.010 1.035 -4.925 -0.993\n", "α1 - α2 -5.014 1.621 -8.095 -2.047\n", "α1 - α3 -7.019 1.617 -9.971 -3.907\n", "α1 - α4 -0.007 1.532 -2.914 2.883\n", "α2 - α3 -2.005 1.438 -4.775 0.659\n", "α2 - α4 5.007 1.334 2.431 7.451\n", "α3 - α4 7.011 1.331 4.489 9.524" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "az.summary(trace, var_names=[\"α\"], filter_vars=\"like\", kind=\"stats\")" ] }, { "cell_type": "code", "execution_count": 10, "id": "74133f95", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Last updated: Wed Mar 22 2023\n", "\n", "Python implementation: CPython\n", "Python version : 3.11.0\n", "IPython version : 8.9.0\n", "\n", "pytensor: 2.10.1\n", "\n", "arviz: 0.14.0\n", "pymc : 5.1.2\n", "numpy: 1.24.2\n", "\n" ] } ], "source": [ "%load_ext watermark\n", "%watermark -n -u -v -iv -p pytensor" ] }, { "cell_type": "code", "execution_count": null, "id": "054298a2", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.0" } }, "nbformat": 4, "nbformat_minor": 5 }