{ "cells": [ { "cell_type": "markdown", "id": "d68537ba", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "(bart_categorical)=\n", "# Categorical regression\n", "\n", ":::{post} May, 2024\n", ":tags: BART, regression\n", ":category: beginner, reference\n", ":author: Pablo Garay, Osvaldo Martin\n", ":::" ] }, { "cell_type": "markdown", "id": "0cf4f392-fdc7-4175-9e72-c8a334abea84", "metadata": {}, "source": [ "In this example, we will model outcomes with more than two categories. \n", ":::{include} ../extra_installs.md\n", ":::" ] }, { "cell_type": "code", "execution_count": 1, "id": "7c087cca", "metadata": {}, "outputs": [], "source": [ "import os\n", "import warnings\n", "\n", "import arviz as az\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import pymc as pm\n", "import pymc_bart as pmb\n", "import seaborn as sns\n", "\n", "warnings.simplefilter(action=\"ignore\", category=FutureWarning)" ] }, { "cell_type": "code", "execution_count": 2, "id": "25cf7b45", "metadata": {}, "outputs": [], "source": [ "# set formats\n", "RANDOM_SEED = 8457\n", "az.style.use(\"arviz-darkgrid\")" ] }, { "cell_type": "markdown", "id": "e73740d8-8e70-48b4-b6f9-eb0c1f7ce72f", "metadata": {}, "source": [ "## Hawks dataset \n", "\n", "Here we will use a dataset that contains information about 3 species of hawks (*CH*=Cooper's, *RT*=Red-tailed, *SS*=Sharp-Shinned). This dataset has information for 908 individuals in total, each one containing 16 variables, in addition to the species. To simplify the example, we will use the following 5 covariables: \n", "- *Wing*: Length (in mm) of primary wing feather from tip to wrist it attaches to. \n", "- *Weight*: Body weight (in gr). \n", "- *Culmen*: Length (in mm) of the upper bill from the tip to where it bumps into the fleshy part of the bird. \n", "- *Hallux*: Length (in mm) of the killing talon. \n", "- *Tail*: Measurement (in mm) related to the length of the tail. \n", "\n", "Also we are going to eliminate the NaNs in the dataset. With these we will predict the \"Species\" of hawks, in other words, these are our dependent variables, the classes we want to predict. " ] }, { "cell_type": "code", "execution_count": 3, "id": "71f3a9bc-979f-44fc-8227-133349e4dfb1", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", " | Wing | \n", "Weight | \n", "Culmen | \n", "Hallux | \n", "Tail | \n", "Species | \n", "
---|---|---|---|---|---|---|
0 | \n", "385.0 | \n", "920.0 | \n", "25.7 | \n", "30.1 | \n", "219 | \n", "RT | \n", "
2 | \n", "381.0 | \n", "990.0 | \n", "26.7 | \n", "31.3 | \n", "235 | \n", "RT | \n", "
3 | \n", "265.0 | \n", "470.0 | \n", "18.7 | \n", "23.5 | \n", "220 | \n", "CH | \n", "
4 | \n", "205.0 | \n", "170.0 | \n", "12.5 | \n", "14.3 | \n", "157 | \n", "SS | \n", "
5 | \n", "412.0 | \n", "1090.0 | \n", "28.5 | \n", "32.2 | \n", "230 | \n", "RT | \n", "
\n", "\n" ], "text/plain": [ "\n" ] }, "metadata": {}, "output_type": "display_data" }, { "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 91 seconds.\n", "Sampling: [y]\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a21c89100700452ea64ed94c28a1c896", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n" ], "text/plain": [] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n" ], "text/plain": [ "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "with model_hawks:\n", " idata = pm.sample(chains=4, compute_convergence_checks=False, random_seed=123)\n", " pm.sample_posterior_predictive(idata, extend_inferencedata=True)" ] }, { "cell_type": "markdown", "id": "fb2e357e-502e-4ac5-9d53-928437bd2a4e", "metadata": {}, "source": [ "## Results \n", "\n", "### Variable Importance \n", "\n", "It may be that some of the input variables are not informative for classifying by species, so in the interest of parsimony and in reducing the computational cost of model estimation, it is useful to quantify the importance of each variable in the dataset. PyMC-BART provides the function {func}`~pymc_bart.plot_variable_importance()`, which generates a plot that shows on his x-axis the number of covariables and on the y-axis the R$^2$ (the square of the Pearson correlation coefficient) between the predictions made for the full model (all variables included) and the restricted models, those with only a subset of the variables. The error bars represent the 94 % HDI from the posterior predictive distribution. " ] }, { "cell_type": "code", "execution_count": 10, "id": "a9d1c616-8c1f-4907-ad5a-adffb290c0c2", "metadata": { "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAAE3CAYAAACq3N6VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA21ElEQVR4nO3deVxVdf7H8fdly1hELFsUlBZBy1xDAVHHpnJCTDI103GNitKmtBqztEUtW61cZsyybMFMxshEGMnRTAyE1HH5lTVaomIuhQvImFzu/f3h496B2C7cC/eAr+fj4eMh53z58r3w4XDe93y/55isVqtVAAAAAGAwHu4eAAAAAABUhrACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJC83D2AmqxatUpbt27V7t279cMPP6ikpERz5szRkCFDatWPxWJRUlKSVqxYoby8PPn6+io6OlqTJ09WSEhIPY0eAAAAQF0ZPqy8+eabys/PV1BQkC677DLl5+fXqZ+nn35aycnJat++vUaPHq1jx44pPT1dmzdv1ieffKLQ0FDXDhwAAACAUww/DWz27Nlav369srOzNWLEiDr1kZ2dreTkZEVEROjTTz/V448/rldeeUULFy7UyZMnNWvWLBePGgAAAICzDH9lJTo62uk+kpOTJUkPP/ywfHx87Nv79eunnj17KjMzU4cPH1br1q2d/loAAAAAXMPwV1ZcYcuWLfL19VX37t0r7OvTp48kKScnp6GHBQAAAKAaTT6sFBcX6/jx4woODpanp2eF/e3atZMk5eXlNfTQAAAAAFSjyYeVwsJCSZK/v3+l+23bbe2qYrFYXDswAAAAANUy/JoVozh16pS7h2B4QUFBOnHihLuHgUaMGoIrUEdwFjUEV6COahYUFFRjmyZ/ZSUgIECSVFRUVOl+23ZbOwAAAADG0OTDiq+vr1q1aqVDhw6ptLS0wn7bWhXb2hUAAAAAxtDkw4ok9ezZU8XFxdq2bVuFfZs2bZIkRURENPSwAAAAAFSjSYWVgoIC7du3TwUFBeW2Dx8+XJL05ptv6ty5c/btGzduVE5OjmJiYtSmTZsGHSsAAACA6hl+gX1ycrK2bt0qSfrhhx/s22zPRenRo4eGDRsmSUpKStKCBQs0adIkPfTQQ/Y+IiMjNWzYMCUnJ2vIkCHq16+fjh8/rrS0NLVo0ULTp09v4FcFAAAAoCaGDytbt25VSkpKuW3btm0rN6XLFlaqM3PmTIWFhWnFihX64IMP5Ovrq1tuuUWTJ09W27ZtXT5uAAAAAM4xWa1Wq7sH0Rhw67macYs+OIsagitQR3AWNQRXoI5qxq2LAQAAADRahBUAAAAAhkRYAQA0CWfOnFHLli1lMpl05swZdw8HAOAChBUAAADARXjjxLUMfzcwo7BYLOJeBNUrLS1VaWmpu4eBRowagjPK1g61BGdQP3AGxyLHWSwWeXhUf+2EsOKgwsJCWSwWdw/D8E6fPu3uIaCRo4ZQV8XFxfb/FxYWcoKAWisuLtYNN9wgSdq1a5d8fX3dPCI0RhyLHOfp6VnjHcEIKw6yWCwymUw1pr8LmZeXlzw9Pd09DDRi1BCcUbZ2PD09qSXUGjUEV6COHGOxWBwKcoSVWvDw8CCsVMPT05PvD5xCDcEZZWuH4zXqghqCK1BHrsV3DwAAAIAhEVYAuF1xcbHatm2r5s2bl5vrCwAALmyEFQAAAACGRFgBAAAAYEiEFQAAAMBFyt7hasuWLdy62EmEFQAAAMAF0tPTNWDAAPvHY8eOVd++fZWenu7GUTVuhBUAQJPAu5lwFjUEZ6SnpysxMVHdunVTVlaWCgsLlZWVpa5duyoxMZHAUkcmq9VqdfcgGoPvv/++ymdAmEwmeXn975E1JSUlVfbjTFuz2ayqflz11VaSvL29HWrbokULnTlzxuX9SucfFmgymSSd/2NisVgM3bZsrdSmbU0PSGpsbT08POwPw6qu7X//+1916dJFpaWl2rNnjy6++GKZzWaH+rVarW5p21C/903pGFGfbb/44gvNnj1bhw4dsm8LDg7Wk08+We5dToljhJHaOnqM+H3b+vi9T09P16xZs6qtIY4RtW8rGeMYUd+/96Wlperbt6+6du2qzz77rNz5osViUXx8vHbs2KGNGzfK09OTY4SHhywWi/z8/HTppZdW2VbioZAOW758eZX7goODdeutt9o//vjjj6s8gF1xxRWKjY21f5ycnKyzZ89W2vbSSy/V7bffbv/4008/VVFRUaVtW7RooSFDhtg//vzzz3Xy5MlK2/r7+2v48OH2j9PS0vTLL79U2rZZs2YaOXKk/eOMjAwdOXKk0rZeXl4aM2aM/eP169eXO+j/3oQJE+z//+qrr7R///4q244ePdp+UNq8ebP27t1bZdu7775bF198saTz74zt2bOnyrbDhg1TQECAJGnr1q3avXt3lW3vuOMOBQUFSZJ27Nihf//731W2HTRokFq1aiVJ+vbbb5Wbm1tl29tuu01XXnmlJGnPnj3Kzs6usu0tt9yikJAQSdK+ffu0adOmKtv2799fV111lSQpLy9PGzZsqLJtnz591L59e0lSfn6+vvjiiyrbRkZG6rrrrpMkHT16tNp3iiIiInTDDTdIkn799VetXr262jF8+eWXkqSTJ08qJSWlyradOnVSz549JUlFRUVKTk6usm2HDh0UHR0tSTp79qw+/vjjKttee+216tu3r6Tzf/g+/PDDKtuGhobqpptusn9cXVuOEefV1zHi22+/VXJysuLi4pScnKxOnTpp9+7dev755zVx4kQNGzbMXrMSxwibxnaM6Nq1q7p37y7J9ceIU6dOKTExscYa4hhxXmM7Rkj1fx6Rk5OjgwcPasWKFRXe2Pbw8NC0adMUHR2tnJwcRUVFcYwoc0yuCWEFANBoWSwWrVu3TnFxceXezYyMjNSqVas0ePBg/etf/1KHDh14ijQqZbFYNHv2bIdqCKjKsWPHJJ0PypWxbbe1g+OYBuYgpoHV3JZpYP9zIVy+daQt08Aatq3RjxH10TY7O1ujRo1SVlaWIiMjK+zPyspSdHS0kpKS7Ps5RhinrRGmgeXk5Ojuu+92qIaio6M5RtSyrXRhTAPLysrSXXfdVWMdffLJJ4qKiuIYwTQw1/P29q4yrFTWtjb9OqrsgcGIbX//WuprDJ6envY/Fk2trYeHh8Pv/jaltiUlJeUObCaTyeHfDSO0lerv974pHSPqo21BQYGkmt/NLCgoqPR7aYTfe44RtW/ryt9l2/QlR2qo7M+JY0TjbVsfv589e/ZUSEiIXnjhhUrXrMyZM0dt27a1T0/kGOE4rokDABqtyy67TJKqXEti225rB/weNQRX8PT01PTp05Wamqr4+PhydwOLj49XamqqnnrqKYdDB/6HsAIAaLTKvpv5+2kSlb2bCfweNQRXue2227Ro0SJt375d0dHRat68uaKjo7Vjxw4tWrRIt912m7uH2CgRVgAAjRbvZsJZ1BBc6bbbbtPatWvtH7///vvauHEjQcUJrFkBADRqtnczZ82aZb9NtSS1bduWdzPhEGoIrlQ22Pbq1Yug6yTCCgCg0bvtttsUExOj66+/XtL5dzP79u3LSQIcRg0BxsQ0MABuV/ZOYFu2bKn2lodAVXg3E86ihgDjIawAcKv09HQNGDDA/vHYsWPVt2/fap96DQAALgyEFQBuk56ersTERHXr1q3cotauXbsqMTGRwAIAwAWOsALALUpLSzV79mzFxcXps88+U2RkpPz9/RUZGanPPvtMcXFxev7555kSBgDABYwF9rXw+/uvo7zS0lK+R3BYdna2Dh48qBUrVlR40q2Hh4emTZum6OhoZWdnKyoqyk2jRGNS9vhjsVg4HqHWqCG4AnXkGEe/L4QVB3l4eMhisfAubzXMZjPfHzjsyJEjkqROnTpVut+2/ciRI9QVHFK2TkpLS6kb1Bo1BFegjhznyE0sCCsOCggIkNVqdfcwDC0oKMjdQ0AjctVVV0mSdu/ercjIyAr7d+/ebW/XvHnzBh0bGqeyf/QCAgLk5+fnxtGgMaKG4ArUkeMCAwNrbENYcdDvp6mgIk9PT27zCIf17t1bbdu21QsvvKDPPvus3O+YxWLRnDlz1K5dO/Xu3Zu6gkPK1gnHI9QFNQRXoI4c58j5NWfgANzC09NTs2bNUmpqquLj48vdDSw+Pl6pqamaOXMmB3kAAC5gXFkB4DaDBg3S0qVLNX36dEVHR9u3t2vXTkuXLtWgQYPcODoAAOBuhBUAbjVo0CD169dPoaGhkqQVK1aof//+XFEBADRKfn5+KigoUFBQkE6cOOHu4TR6TAOD086cOaOWLVvKZDLpzJkz7h4OGqGywSQqKoqgAgAAJDWSKys7d+7U/PnztX37dpnNZoWFhWncuHGKjY11uI+jR4/q7bff1tdff63Dhw/L19dX7dq101133aVBgwZxcgQAAAAYjOHDSnZ2thISEuTj46OBAwfKz89PGRkZmjx5so4cOaIJEybU2MfBgwc1bNgwnTx5UjExMerfv7+Kior0r3/9S1OnTtWWLVs0Z86cBng1AAAAABxl6LBiNps1Y8YMmUwmJSUlqWPHjpKkiRMnaujQoZo7d64GDBigNm3aVNvPkiVLdOLECT355JMaO3asffujjz6qwYMH69NPP9WkSZNq7AcAAABAwzH0mpXs7GwdOHBAcXFx9qAinX/ATmJiokpKSpSSklJjPwcPHpQk9evXr9z25s2bq3v37pLEAigAaORsi1qtVisPYQOAJsLQYSUnJ0eSFBMTU2GfbVtubm6N/YSFhUmSNm7cWG776dOntX37drVq1UrXXnuts8MFAACNGIEXMB5DTwPbv3+/pPPPXPi9Vq1aydfXV3l5eTX2c88992j9+vWaM2eONm3apPDwcPualWbNmmnBggVq1qyZq4cPAAAAwAmGDitFRUWSzk/7qoy/v78KCwtr7OfSSy/VJ598oscff1xfffWVNm3aJElq1qyZRowYoQ4dOtTYR2BgoDw8DH0hym18fHzs/w8KCuLdKNQaNQRXCwoKcvcQ0MhRQ3AF6sh5hg4rrpKXl6fExET5+vraF+oXFhbq888/1xtvvKHMzEwlJSVVe/viU6dONeCIG5eyz1Y5ceKEzp0758bRoDGihuBKPIgNzqKG4ArUUc0cCXOGDiv+/v6SVOXVk6KiIgUGBtbYzxNPPKHDhw9r3bp1atWqlaTz81Lvu+8+/fLLL3r//fe1Zs0a3X777a4bPAAAAACnGHpeU2hoqCRVui7l+PHjKi4urnQ9S1lFRUXatm2brrnmGntQKatXr16SpO+++875AQMAAABwGUOHlYiICElSZmZmhX22bbY2VSkpKZFU9a2JCwoKJJWfMw8AAADA/QwdVqKiohQSEqLU1NRyVz4KCwu1aNEieXt7Kz4+3r792LFj2rdvX7lpY0FBQbrqqqt0+PBhJScnl+v/9OnTevfddyX97woLAAAAAGMw9JoVLy8vzZ49WwkJCRo1apQGDhwoPz8/ZWRkKD8/X1OnTlVwcLC9/dy5c5WSkqI5c+ZoyJAh9u3Tpk3Tgw8+qOnTp2vNmjXq2LGjTp8+rfXr16ugoEADBgxQdHS0O14iAAAAgCoYOqxIUmRkpJYtW6Z58+YpLS1NZrNZYWFheuyxxxQbG+tQH/369dPHH3+sJUuWaOvWrcrNzZWPj4+uueYaTZw4UXfffXc9vwoAAAAAtWWyWq1Wdw+iMeDWc1U7c+aMQkJCJEkHDx7kGRmoNWoIrsTtQuEsagiuQB3VzJFbFxt6zQoAAACACxdhBQAAAIAhEVYAAAAAGBJhBQAAAIAhEVYAAAAAGBJhBQAAAIAhEVYAAAAAGBJhBQAAAIAhEVYAuJ2fn58KCgpktVp5ICQAALAjrMBppaWl9v9nZWWV+xgAAACoK8IKnLJ69Wr16dPH/vHw4cPVo0cPrV692o2jAgAAQFNAWEGdrV69WuPGjVPXrl2VlZWlwsJCZWVlqUuXLho3bhyBBQAAAE4xWa1Wq7sH0RicOHHC3UMwlNLSUvXo0UNdunTRZ599Jg+P/+Vei8Wi+Ph47dy5U9988408PT3dOFI0JkFBQfyuwWnUEZxFDcEVqKOaBQUF1diGKyuok6ysLB04cEBPPvlkuaAiSR4eHpo2bZry8vKUlZXlphECAACgsSOsoE6OHj0qSerUqVOl+23bbe0AAACA2iKsoE4uv/xySdLu3bsr3W/bbmsHAAAA1BZhBXUSFRWltm3b6oUXXpDFYim3z2KxaM6cOWrXrp2ioqLcNEIAAAA0doQV1Imnp6dmzZql1NRUxcfHl7sbWHx8vFJTUzVz5kwW1wMAAKDOvNw9ADRegwYN0tKlSzV9+nRFR0fbt7dr105Lly7VoEGD3Dg6AAAANHaEFThl0KBB6tevn0JDQyVJK1asUP/+/bmiAgAAAKcxDQxOKxtMoqKiCCoAAABwCcIKAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwJC93D8ARO3fu1Pz587V9+3aZzWaFhYVp3Lhxio2NrVU/v/76q9566y19+eWX+vnnn+Xr66vQ0FANHjxYI0eOrKfRAwAAAKgLw4eV7OxsJSQkyMfHRwMHDpSfn58yMjI0efJkHTlyRBMmTHCon++++04TJkzQ6dOn1a9fPw0YMEDFxcXat2+fNmzYQFgBAAAADMbQYcVsNmvGjBkymUxKSkpSx44dJUkTJ07U0KFDNXfuXA0YMEBt2rSptp+ioiI9+OCDkqSVK1eqQ4cOFb4OAAAAAGMx9JqV7OxsHThwQHFxcfagIkkBAQFKTExUSUmJUlJSauxn2bJlOnz4sB599NEKQUWSvLwMndkAAACAC5Khz9JzcnIkSTExMRX22bbl5ubW2E9aWppMJpMGDBigH3/8UZs3b9bZs2d19dVXq0+fPvLx8XHtwAEAAAA4zdBhZf/+/ZKkdu3aVdjXqlUr+fr6Ki8vr9o+zp07px9++EEtW7bUhx9+qPnz58tisdj3h4SEaOHChQoPD3fp2AEAAAA4x9BhpaioSNL5aV+V8ff3V2FhYbV9nDp1SqWlpTp58qT+9re/6fHHH9fgwYNlNpu1fPly/f3vf9cDDzyg9PR0XXTRRVX2ExgYKA8PQ8+ac5uyV6aCgoLk5+fnxtGgsQsKCnL3ENAEUEdwFjUEV6COnGfosOIKtqsopaWlGjVqVLm7hz388MP66aeflJ6ern/+858aPHhwlf2cOnWq3sfaWJ05c8b+/xMnTujcuXNuHA0as6CgIJ04ccLdw0AjRx3BWdQQXIE6qpkjYc7Qlwr8/f0lqcqrJ0VFRVVedbEpu/+mm26qsN+2bffu3XUdJgAAAIB6YOiwEhoaKkmVrks5fvy4iouLK13PUpavr68uv/xySVLz5s0r7Ldt++2335wcLQAAAABXMnRYiYiIkCRlZmZW2GfbZmtTncjISEnS3r17K+yzbavpWS2omp+fnwoKCmS1WlmvAgAAAJcxdFiJiopSSEiIUlNT9d1339m3FxYWatGiRfL29lZ8fLx9+7Fjx7Rv374K08ZGjBghSXr77bd1+vRp+/bjx4/rgw8+kIeHh2699db6fTEAAAAAasXQC+y9vLw0e/ZsJSQkaNSoURo4cKD8/PyUkZGh/Px8TZ06VcHBwfb2c+fOVUpKiubMmaMhQ4bYt3fv3l3jx4/Xe++9p9tvv139+/eX2WzWv/71L/3666+aMmWKrrrqKne8RAAAAABVqHVY+f3zSEwmk/z8/HTNNdcoNjZWo0aNkre3d7k2aWlpWrp0qQ4dOiSz2az27dvr0UcfVffu3Wv8epGRkVq2bJnmzZuntLQ0mc1mhYWF6bHHHlNsbKzD437iiScUFhampKQkpaSkyGQyqWPHjnruued0yy23ONwPAAAAgIZhslqt1tp8gi2s3HHHHZLO3xI4Pz9f27dvl8ViUVRUlN555x15ef0vB6WlpSkoKEg33nijSktL9dxzz2nt2rVKT0+3L343Om49VzNu0QdnUUNwBeoIzqKG4ArUUc0cuXVxnaeBvfjii+U+3rFjh0aPHq2srCytWbOm3DNLyl4B8fb21tixY/Xpp59q165djSasAAAAAGhYLltg36VLF/vVlsru3lXW2rVr5eXlVWFKGQAAAADYuPRuYO3bt5ckFRQUVNnmm2++0dtvv60xY8YoJCTElV8eAAAAQBPi0rBy5swZSVLLli0r3b9v3z5NnDhRERERevTRR135pQEAAAA0MS4NK5s2bZIk9enTp8K+AwcOaNy4cWrfvr0WLlxYbgE+AAAAAPye04nBYrHo0KFDWrJkiXJzc/XHP/6xwi2FDx48qDFjxujaa6/VwoUL5evr6+yXBQAAANDE1TmsVLY4fvjw4Zo5c6ZMJlO57X/729/0888/68yZM+rfv799e0JCgu699966DgEAAABAE+b0c1Z+++037dmzRz/++KMkVXh6fFPBfbJrxv3E4SxqCK5AHcFZ1BBcgTqqWYM+Z+Wdd97RK6+8opkzZ6pXr15q06ZNXbsGAAAAANctsE9ISFBMTIz++9//asGCBa7qFgAAAMAFyqV3A7Pdjvjzzz9Xfn6+K7sGAAAAcIFxaVi57rrrdPPNN8tsNuudd95xZdcAAAAALjAuDSuS9NBDD8lkMmnlypU6fvy4q7sHAAAAcIFweVjp0KGDbrnlFv3222967733XN09AAAAgAtEre8G9v3339fYZv78+XUaDAAAAADYuPzKCgAAAAC4AmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCE1irCyc+dO3XvvvbrxxhvVtWtXDR8+XGlpaXXu79SpU+rTp4/Cw8N1zz33uHCkAAAAAFzFy90DqEl2drYSEhLk4+OjgQMHys/PTxkZGZo8ebKOHDmiCRMm1LrPmTNnqqioqB5GCwAAAMBVDH1lxWw2a8aMGTKZTEpKStKsWbP0xBNPaNWqVQoNDdXcuXOVn59fqz7Xrl2r1NRUPfbYY/U0agAAAACuYOiwkp2drQMHDiguLk4dO3a0bw8ICFBiYqJKSkqUkpLicH8FBQV69tlnNXjwYPXr168+hgwAAADARQwdVnJyciRJMTExFfbZtuXm5jrc3zPPPCNPT0899dRTrhkgAAAAgHpj6DUr+/fvlyS1a9euwr5WrVrJ19dXeXl5DvW1atUqZWRkaOHChQoMDFRhYaErhwoAAADAxQwdVmyL4AMCAird7+/v71DoOHr0qJ5//nnFxcXp5ptvrtNYAgMD5eFh6AtRhhAUFOTuIaCRo4bgCtQRnEUNwRWoI+cZOqy4yvTp0+Xl5eXU9K9Tp065cERNU1BQkE6cOOHuYaARo4bgCtQRnEUNwRWoo5o5EuYMHVb8/f0lqcqrJ0VFRQoMDKy2j5SUFH311Vd688031bJlS5ePEQAAAED9MHRYCQ0NlSTl5eWpU6dO5fYdP35cxcXF6ty5c7V9fPvtt5Kkhx9+uNL9mZmZCg8PV4cOHbRq1SrnBw0AAADAJQwdViIiIvTWW28pMzNTAwcOLLcvMzPT3qY63bp1U3FxcYXtxcXFSktL0xVXXKGYmBhdeeWVrhs4AAAAAKcZOqxERUUpJCREqampGjNmjP1ZK4WFhVq0aJG8vb0VHx9vb3/s2DEVFhbqsssusy/Kj42NVWxsbIW+Dx06pLS0NF177bV6/vnnG+T1AAAAAHCcoW9v5eXlpdmzZ8tqtWrUqFGaMWOGXnzxRQ0ePFj79+/XlClTFBwcbG8/d+5cxcbG6osvvnDjqAEAAAC4gqGvrEhSZGSkli1bpnnz5iktLU1ms1lhYWF67LHHKr1iAgAAAKBpMFmtVqu7B9EYcOu5mnGLPjiLGoIrUEdwFjUEV6COaubIrYsNPQ0MAAAAwIWLsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAyJsAIAAADAkAgrAAAAAAzJy90DcMTOnTs1f/58bd++XWazWWFhYRo3bpxiY2Nr/Fyr1aqvvvpK69ev17Zt23T48GGZzWa1a9dOsbGxGj9+vC666KIGeBUAAAAAasPwYSU7O1sJCQny8fHRwIED5efnp4yMDE2ePFlHjhzRhAkTqv38c+fO6b777pOPj4969uypmJgYnTt3TpmZmXr99de1bt06ffjhh7r44osb6BUBAAAAcIShw4rZbNaMGTNkMpmUlJSkjh07SpImTpyooUOHau7cuRowYIDatGlTZR8eHh565JFHNHLkSAUGBtq3l5SU6KGHHtKGDRuUlJSkhISEen89AAAAABxn6DUr2dnZOnDggOLi4uxBRZICAgKUmJiokpISpaSkVNuHt7e3HnjggXJBxbb9/vvvlyTl5ua6fvAAAAAAnGLosJKTkyNJiomJqbDPts2ZoOHldf7CkqenZ537AAAAAFA/DB1W9u/fL0lq165dhX2tWrWSr6+v8vLy6tz/ypUrJUm9e/eucx8AAAAA6oeh16wUFRVJOj/tqzL+/v4qLCysU98bN27UJ598omuuuUbDhg2rsX1gYKA8PAyd7QwhKCjI3UNAI0cNwRWoIziLGoIrUEfOM3RYqS87d+7U5MmTFRAQoDfffFM+Pj41fs6pU6caYGSNW1BQkE6cOOHuYaARo4bgCtQRnEUNwRWoo5o5EuYMfanA399fkqq8elJUVFTlVZeq7Nq1S/fcc488PDz0zjvvqH379k6PEwAAAIDrGTqshIaGSlKl61KOHz+u4uLiStezVGXXrl2aMGGCLBaLlixZos6dO7tqqAAAAABczNBhJSIiQpKUmZlZYZ9tm61NTWxBpbS0VO+88466dOniuoECAAAAcDlDh5WoqCiFhIQoNTVV3333nX17YWGhFi1aJG9vb8XHx9u3Hzt2TPv27aswbWz37t2aMGGCzGaz3n77bXXr1q2hXgIAAACAOjL0AnsvLy/Nnj1bCQkJGjVqlAYOHCg/Pz9lZGQoPz9fU6dOVXBwsL393LlzlZKSojlz5mjIkCGSpJMnT2rChAk6ffq0+vTpo6+//lpff/11ua8TEBCgcePGNeRLAwAAAFADQ4cVSYqMjNSyZcs0b948paWlyWw2KywsTI899phiY2Nr/PyioiL7nbw2bdqkTZs2VWjTpk0bwgoAAABgMCar1Wp19yAaA249VzNu0QdnUUNwBeoIzqKG4ArUUc0a/a2LAQAAAFy4CCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQvNw9AEfs3LlT8+fP1/bt22U2mxUWFqZx48YpNjbW4T7OnTunxYsX6/PPP9fPP/+swMBA9e/fX4888oguueSSehw9AAAAgLowfFjJzs5WQkKCfHx8NHDgQPn5+SkjI0OTJ0/WkSNHNGHChBr7sFgseuCBB5SZmamuXbvq1ltvVV5enpKTk5WVlaUVK1aoZcuWDfBqAAAAADjK0GHFbDZrxowZMplMSkpKUseOHSVJEydO1NChQzV37lwNGDBAbdq0qbaflJQUZWZmKi4uTq+++qpMJpMk6eOPP9azzz6rN954QzNnzqz31wMAAADAcYZes5Kdna0DBw4oLi7OHlQkKSAgQImJiSopKVFKSkqN/SQnJ0uSpkyZYg8qkjRixAiFhIRo9erVOnv2rOtfAAAAAIA6M3RYycnJkSTFxMRU2GfblpubW20fv/32m3bs2KGrrrqqwhUYk8mk6OhoFRcXa/fu3S4aNQAAAABXMHRY2b9/vySpXbt2Ffa1atVKvr6+ysvLq7aPAwcOyGKxKDQ0tNL9tu22rwUAAADAGAy9ZqWoqEjS+WlflfH391dhYWG1fdj2+/v7V9lH2a9VlaCgoGr34zy+T3AWNQRXoI7gLGoIrkAdOc/QV1YAAAAAXLgMHVZsVz2qunpSVFRU5VUXG9v+qq6c2LZXdeUFAAAAgHsYOqzY1pNUti7l+PHjKi4urnQ9S1khISHy8PCock2KbXtVa1oAAAAAuIehw0pERIQkKTMzs8I+2zZbm6o0a9ZMnTt31k8//aT8/Pxy+6xWq77++mv5+vqqU6dOLho1AAAAAFcwdFiJiopSSEiIUlNT9d1339m3FxYWatGiRfL29lZ8fLx9+7Fjx7Rv374K08aGDx8uSZo7d66sVqt9+/Lly3Xw4EENGjRIzZo1q98XAwAAAKBWDB1WvLy8NHv2bFmtVo0aNUozZszQiy++qMGDB2v//v2aMmWKgoOD7e3nzp2r2NhYffHFF+X6ueOOOxQTE6PU1FSNGDFCr776qv7yl7/oueeeU3BwsB555JEGfmVNx+jRoxUeHl5u25YtWxQeHq758+e7aVRo7ObPn6/w8HBt2bLF3UNBE0NtwVVc9beOmkRNLvRzLUOHFUmKjIzUsmXL1L17d6Wlpenjjz/WJZdcotdff10TJkxwqA8PDw/9/e9/10MPPaSCggItXbpU27Zt09ChQ/XJJ5+oZcuW9fwqGt6hQ4cUHh6ue+65p8o2tkJ/+umnG3BkaAqsVqsyMjI0adIk9e3bV506dVK3bt10++2364UXXtDevXvdPUQ0EYsXL1Z4eLgWL15c6f7BgwdXexx78MEHFR4erqysrPocJtzkQqyPC+kktTGx/Vwc/Td69Gh3D7nRMPRzVmw6d+6sd955p8Z2L774ol588cVK9/n4+GjSpEmaNGmSq4cHXFBOnjyphx9+WNnZ2WrevLmio6MVEhKikpIS7d27V8uWLdOHH36opUuXqlevXu4eLho5Ww1t2bJF9913X7l9J06c0Pfffy+TyaScnJwKn2uxWPTNN9/Ix8dH3bt3lySNGjVKsbGxat26df0PHvXO1fVRG507d1ZaWhrP0YAkqU2bNhXOMU+fPq0PPvhAbdq00R133FGhvaNeeukl/fe//3XJOBujRhFWADjviSeeUEpKir7//vs692E2mzVp0iTl5ubq9ttv1zPPPFPhtt/Hjh3T66+/XuMDW3FhcLburr/+evn5+Wnbtm0ym83y8vrfn63c3FxZrVbdeuutysjI0LFjx3TZZZfZ9+/Zs0enTp1Sz549ddFFF0mSWrZs2SSvpjdWRquP2rj44ot1zTXX1GncMCZn6jE4OFgPPfRQuW2HDh2yh5Xf76uNC/3NFcNPA0PD2r17t2bOnKm4uDj16NFDnTt31qBBg7R48WKVlJQ41Xd1lz1vuukm3XTTTfaP8/Ly1K1bN/Xt21cnTpwo17a6fahfq1atUm5uriIiIvTSSy9V+nyiyy67THPmzFHfvn3t22rzs6+KbWrjE088oX379un+++/XjTfeqIiICE2ZMkUFBQWSpO3bt2vs2LHq3r27IiIi9NRTT6m4uLjSPnNzc5WYmKhevXqpU6dOuvXWW/X6669XeAer7LSLXbt2afz48erWrZt69OihiRMn6tChQzWOH3Xj5eWlG2+8UcXFxdq1a1e5fTk5OWrWrJnuvfdeSaow59/2bnrZK3yVrQ8oW1t5eXmaOHGiIiIi1LVrV40bN0579uypdGw5OTkaNWqUunbtql69eumRRx7Rzz//XOn8ctQPV9eHJK1bt05jx45VRESEbrjhBsXFxWnJkiUqLS0t16666Vh1rY3Vq1dr8ODB6ty5s2JiYjR79mydPXvWvn/+/PkaM2aMJGnBggXlphVxHGo8anuudaEfUwgrKGfFihX64osvFBYWprvuuktDhw6V1WrVa6+9pilTpjTYONq1a6cZM2bo6NGjmj59un17SUmJpkyZorNnz+rll1/m8nsD+8c//iFJeuCBB+ThUf3hw8fHp17GcOjQIY0YMULnzp3TsGHD1KFDB61Zs0YTJ07UN998o3HjxsnX11d33XWXQkJC9I9//EOzZs2q0M+yZcs0evRobdu2TX/4wx80evRoXX755Vq0aJHGjx+vc+fOVficXbt26c9//rO8vb01YsQIderUSevWrdP48eP122+/1cvrRfmpPmVt2bJFXbp0UadOnRQYGFhhf3Z2drnPr0l+fr6GDx+uU6dO6c4771R0dLSysrI0ZswY/fLLL+XaZmZmavz48dq5c6cGDBig4cOH6/Dhwxo5cqROnz5d15eKOnBlfbz22muaOHGifvrpJ91yyy0aOXKkLrroIr388suaPHmyQ+Opa20kJSXp6aef1rXXXqu7775bzZs314cffqinnnrK3qZnz5726UQ9e/a0T2+fNGmSmjdv7tD44H5GOddqLJgG1sQdOHCgykV4v3/ujCQlJibqmWeekaenp32b1WrVU089pZUrV2rr1q3q0aNHvY23rCFDhigzM1Nr1qzRsmXLNHLkSL3++uvavXu37r//fkVGRjbIOHCe2WzWrl277O9kuktubq6efPJJjR07VtL5+rz//vu1ceNGPfDAA3rttdd08803Szofbu+88059/vnnevTRR3XppZdKkvbu3avnn39e4eHhWrp0abnQu3jxYr322mv66KOPKtzEY+PGjXr99dcVGxtr3/bXv/5Vq1at0rp16zRw4MD6fvkXJNvJZE5OjhITEyVJBQUF+s9//qNJkybJw8NDPXr0KHcyarFYtHXrVjVr1kxdunRx6Ovk5OTo0UcfLbf24Y033tDf//53ffrpp/btpaWlevrpp1VaWqr333+/3O/D1KlT9dlnnzn7klELrqqPzZs3a/HixYqJidH8+fPl6+sr6fwx5tlnn9Xy5cu1du1aDRgwoMqxOFMbX3/9tVauXKmrr75akjR58mQNHjxYaWlp+utf/6rLL7/c/lpTUlLUs2dPp6YWwX2MdK7VGHBlpYk7cOCAFixYUOm/lJSUCu1bt25d7pdHkkwmk0aNGiVJDX7HlOeee05t2rTRSy+9pA8//FDvvvuuOnfurL/85S8NOg6cX1hfUlKioKCgOs3vdpW2bdvap0FI5+vTFh46duxoDyqS5O3trQEDBshsNpe7Q9ny5ctlNps1Y8aMClfnEhIS1LJlS6Wmplb42hEREeWCiiTdeeedklRhCgpc57rrrlPz5s21bds2+xWvnJwcWa1W+8lbz549deDAAf3888+SpO+++06nT59Wt27dHL7KFxwcrISEhHLbhg4dKqn8z3fr1q3Kz89X//79KwT3Rx55pMIxFPXLVfXx0UcfSZJmzZplDyrS+WPMY489JpPJpDVr1lQ7FmdqY8yYMfagIp1/qHVcXJwsFov+7//+z9FvBxoBo51rGR1XVpq4mJgYLVmypNJ9W7ZsKXfSJ0nnzp1TUlKS1qxZox9//FHFxcXlHqR57Nixeh3v7wUEBOjVV1/Vn//8Z82ePVt+fn567bXXyi2iREU33XRTpVfOJFU673XOnDkaMmRIfQ/LJcLDw2Uymcptsy2a7dixY4X2tn1la3fHjh2SpE2bNlX6R8HLy0s//fRThe3XX399hW1XXHGFJDH1R/VXd7Z3xjds2KCdO3fqxhtv1JYtW3TRRRfZ3xWPiIiQdP64Fh8fX+V6hOp07NixwvTGyn6+tjUslb3zeeWVV+rKK69k/UAljF4fO3bskK+vr1auXFnp12nWrJl+/PHHasfiTG1wfGlY7vw7abRzLaPjjA/l/OUvf9GGDRsUGhqq2NhYXXLJJfLy8rLffq+yefz17frrr1fr1q118OBB9e3bV23btm3wMTQ2Y8aMqXA3rnXr1mnPnj2V3r67spP832vRooW8vb118uRJnTt3rt7WpNSkskX9tneoqttnNpvt206dOiVJWrRokcu+tsViqVVfTVF91J1NZGSkNmzYoJycHN14443KyclRly5d7HXYsWNHBQQE2E9GbVN+ajNdtLKfr+2NkbI/36KiIknSJZdcUmk/l156KWGlEkavj1OnTslsNmvBggVVfp2qbtZh40xtcHxpWPVZjzUx4rmWkRFWYLdz505t2LBBMTExWrx4cblLlP/+97/1wQcfONW/yWQqd8JYVmFhoQICAird9/LLL+vgwYNq0aKF0tPTdccdd6hfv35OjaWpGzduXIVt+fn52rNnT53nOHt5eemGG27Qtm3blJubq969ezv8uXX92dcX20nB1q1bKz1BQN3UR93ZlF1Efdddd2nv3r3l+vT09LSvS7A9P8PX11c33HCDU1+3Mraa+fXXXyvd//vF+DjP6PVh+7k68yR5aqPxqM96rE59n2s1RaxZgd3BgwclSX/4wx8qzKX85ptvnO4/MDBQR48erbD90KFDVV7i/vLLL/XRRx+pZ8+eWrlypQIDAzVt2jQO+G5im7+/aNGicpesK1P2naG6/OzrU+fOnSX9bzoYjC88PFwtWrTQv//9b23atEnS+XUIZUVERCg/P18ZGRkqLCxUjx496mXKaIcOHSRJ27Ztq7DvyJEj9nURaDiuqI/OnTvr5MmT2r9/f53H0RC1Yfv7/PtbKaNxqO9zraaIsAI720OHtm7dWm77f/7zHy1evNjp/jt16qT8/PxyTxI+d+6cXnzxxUrbHz9+XNOmTVNgYKBeeeUVBQcHa+bMmfr11181derUGk+W4XqDBw+2T7GYNm2afcpDWb/88oumT5+ur776yr6ttj/7+jZy5Eh5eXlp1qxZOnz4cIX9p0+f1rfffuuGkaEqHh4eioiI0NmzZ7VkyRJddNFF6tq1a7k2tpPThQsXSqrdepXa6NGjh1q3bq0NGzZo+/bt5fa9+eabnES6gSvqw/YsqCeffLLSZ3gdP35c+/btq3YcDVEbgYGBks6HHzQ+9X2u1RQxDQx2nTt3VufOnZWenq7jx4+rS5cu+vnnn7V+/Xr169dPa9eudar/8ePHa/Pmzbrvvvs0cOBAXXzxxdq8ebOaN2+uVq1alWtrtVo1depUFRQUaN68efZFhn/60580dOhQ/eMf/9B7771X4dayqF9eXl5auHChHn74YaWkpGj9+vXq3bu3goODVVJSor179yonJ0dms1m33367/fNq87NvCGFhYXrmmWf07LPP6k9/+pP69eunkJAQnTlzRocOHVJOTo7uuOMOzZw5s8HHhqr16tVLX3zxhX744Qf17Nmzwrqp6667Tr6+vvrhhx/s7euDp6ennn32WT344IMaO3asYmNj1apVK+Xm5uro0aPq0KFDnZ/Ijrpztj769u2rBx98UH/729906623qk+fPmrdurVOnjypvLw8bd26VY888ki1T61viNq4+uqrddlll2nNmjXy8fHR5ZdfLpPJpNGjRzf4lFrUXn2fazVFXFmBnaenp9566y3deeedOnDggD766CPt3btXf/3rX/X444873X9MTIzeeOMNhYSEaNWqVfrnP/+p3r176913363wR+Xdd9/V5s2bNWzYsAr3tH/qqacUGhqquXPn8u63G7Ro0UJLly7VvHnzFBERoW+++Ubvvfeeli9frqNHj2r48OFavXp1uSkYtfnZN5Thw4dr+fLluvnmm+3zhNeuXasTJ05o3Lhx9ue4wDjKnlxWFkS8vLzUvXt3SefXDlR2dyVX6devn5YsWaJOnTopPT1dK1as0OWXX65ly5bJYrGwFsoNXFEfDz/8sN577z316NFDWVlZWrp0qb788kuVlJRo0qRJGjRoUI3jqO/a8PT01IIFC9S1a1elpqZq3rx5evPNN+03DoGx1fe5VlNksjKXBgAAlygqKlLv3r0VFham5ORkdw8HBkJtAHXDlRUAAGqpuLi4wpqt0tJSvfzyyzp79my5h5PiwkJtAK7FmhUAAGopLy9PI0eOVExMjIKDg3XmzBlt3bpVe/fuVfv27e2LtXHhoTYA12IaGAAAtVRQUKBXXnlFOTk5+vXXX2U2m9W6dWv98Y9/1AMPPKDmzZu7e4hwE2oDcC3CCgAAAABDYs0KAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwJMIKAAAAAEMirAAAAAAwpP8Ho157Tj9o36UAAAAASUVORK5CYII=", "text/plain": [ "
<xarray.DataArray 'y' ()> Size: 8B\n", "array(96.56689113)
<xarray.DataArray 'y' ()> Size: 8B\n", "array(96.56689113)
\n", "\n" ], "text/plain": [ "\n" ] }, "metadata": {}, "output_type": "display_data" }, { "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 262 seconds.\n", "Sampling: [y]\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f07ff8fa27e647099008bfbf163d8e4f", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n" ], "text/plain": [] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n" ], "text/plain": [ "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "with pm.Model(coords=coords) as model_t:\n", " μ_t = pmb.BART(\"μ\", x_0, y_0, m=50, separate_trees=True, dims=[\"species\", \"n_obs\"])\n", " θ_t = pm.Deterministic(\"θ\", pm.math.softmax(μ_t, axis=0))\n", " y_t = pm.Categorical(\"y\", p=θ_t.T, observed=y_0)\n", " idata_t = pm.sample(chains=4, compute_convergence_checks=False, random_seed=123)\n", " pm.sample_posterior_predictive(idata_t, extend_inferencedata=True)" ] }, { "cell_type": "markdown", "id": "60dc23a9-2351-4502-824a-944e0f454c4c", "metadata": {}, "source": [ "Now we are going to reproduce the same analyses as before. " ] }, { "cell_type": "code", "execution_count": 16, "id": "a05a3d39-307a-4c08-93ec-3a0503ea6c25", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAAE3CAYAAACq3N6VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3RUlEQVR4nO3deVzU1eL/8fewaSwillYKQosgRe4oGuq1W3lFTDI10+saFaUt2mJuLWrZqrmVWZYtmmmGJsKVupqJiZBa6resqyUo5pKogFyTYfj94W/mSuzMwHzA1/Px6PGQcw6HM3Q483nP53w+H1NRUVGRAAAAAMBgXJw9AAAAAAAoDWEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCERVgAAAAAYEmEFAAAAgCG5OXsAFVm7dq127NihvXv36pdfflFBQYFmzZqlAQMGVKkfi8WiZcuWaeXKlcrIyJCnp6e6deum8ePHKyAgoIZGDwAAAKC6DB9W5s6dq6ysLPn5+alZs2bKysqqVj/PPPOMVq1apVatWmn48OE6fvy4kpKStHXrVn366acKCgpy7MABAAAA2MXw28BmzpypjRs3KjU1VUOGDKlWH6mpqVq1apXCw8P1+eef68knn9Srr76qhQsX6vTp05oxY4aDRw0AAADAXoY/s9KtWze7+1i1apUk6dFHH5WHh4etvGfPnurcubNSUlJ05MgRNW/e3O6fBQAAAMAxDH9mxRG2b98uT09PdejQoURd9+7dJUlpaWm1PSwAAAAA5aj3YSU/P18nTpyQv7+/XF1dS9QHBgZKkjIyMmp7aAAAAADKUe/DSm5uriTJ29u71HprubVdWSwWi2MHBgAAAKBchr9mxSjOnDnj7CEYnp+fn06dOuXsYaAOYw7BEZhHsBdzCI7APKqYn59fhW3q/ZkVHx8fSVJeXl6p9dZyazsAAAAAxlDvw4qnp6eaNm2qw4cPq7CwsES99VoV67UrAAAAAIyh3ocVSercubPy8/O1c+fOEnVbtmyRJIWHh9f2sAAAAACUo16FlezsbB04cEDZ2dnFygcPHixJmjt3rs6fP28r37x5s9LS0hQZGakWLVrU6lgBAAAAlM/wF9ivWrVKO3bskCT98ssvtjLrc1E6duyoQYMGSZKWLVumBQsWaNy4cXr44YdtfURERGjQoEFatWqVBgwYoJ49e+rEiRNKTExU48aNNXXq1ArHYbFYVFRU5OiXV68UFhaWutUOqCzmEByBeQR7MYfgCMyjilksFrm4lH/uxPBhZceOHYqPjy9WtnPnzmJbuqxhpTzTp09XcHCwVq5cqQ8//FCenp667bbbNH78eLVs2bLC78/NzeX2xZWQk5Pj7CGgjmMOwRGYR6iO/Px83XTTTZKkPXv2yNPT08kjQl3HWlQ+V1fXCu8IZiridEGl/PrrrzKZTBWmv0uZr68vt3iGXZhDsEd+fr5uuOEGSdKPP/7IgSaqjDkER+I9rXwWi0U+Pj664oorym1n+DMrRuLi4kJYKYerqyu/H9iFOQR7XDx3WK9RHcwhOEJ+fr5at24tSdq3bx+h1078FQIAAAAwJMIKAAAAAEMirAAAAEjF7ty0fft27uQEGABhBYDT5efnq2XLlmrUqJHy8/OdPRwAl6CkpCT17t3b9vXIkSPVo0cPJSUlOXFUqIsIvY5FWAEAAJe0pKQkxcXFqX379tq2bZtyc3O1bds2tWvXTnFxcQQWVBqh1/EIKwCAeoFPM1EdhYWFmjlzpqKjo7VmzRpFRETI29tbERERWrNmjaKjo/XCCy8wn1AhQm/NIKwAAOo8Ps1EdaWlpenQoUOaPHlyiVsVu7i4aNKkScrMzFRaWpqTRoi6gNBbc3jOSiUVFBTIYrGUes91k8kkNze3Ym3LYk9bs9mssp7hWVNtJcnd3b1Sbf/6WhzVryS5ubnJZDJJurAgWCwWQ7e9+HkhVWlrsVjKXcjqWlsXFxe5urpW2NZsNsvV1dVWX1RUJLPZXKl+ndW2tv7u69MaUVNtN2zYoLFjxyo6OlqrVq1SWFiY9u7dqxdffFFxcXFauHBhsSDDGmGctpVdI/7a1pF/y0ePHpUkhYWFlVpvLT9y5IgKCwtZI6rYVnL+GiHV/N+9NfSuXLmyzNDbrVs3paWlqWvXrqwRVXiGEWGlklasWFFmnb+/v26//Xbb15988kmZC9hVV12lqKgo29erVq3SuXPnSm17xRVX6I477rB9/fnnnysvL6/Uto0bN9aAAQNsX3/xxRc6ffp0qW29vb01ePBg29eJiYn6448/Sm3bsGFDDR061PZ1cnKybWH/Kzc3N40YMcL29caNG3X48OFS20rSmDFjbP/+5ptvdPDgwTLbDh8+3LYobd26Vfv37y+z7T333KPLLrtM0oWtIPv27Suz7aBBg+Tj4yNJ2rFjh/bu3Vtm2zvvvFN+fn6SpB9++EHff/99mW379eunpk2bSrrwFOT09PQy2/bp00dXX321pAsPj0pNTS2z7W233aaAgABJ0oEDB7Rly5Yy2/bq1UvXXHONJCkjI0ObNm0qs2337t3VqlUrSVJWVpa+/PLLMttGRETYnvB87Nixcj+5Dg8P10033SRJOnnypNatW1fuGL7++mtJ0unTpxUfH19m27CwMHXu3FmSlJeXp1WrVpXZtnXr1urWrZsk6dy5c/rkk0/KbHv99derR48eki688X300Udltg0KCtItt9xi+7q8tqwRF9TEGmGxWLRgwQLbp5nWN0Drp5n9+/fXtGnTdPToUVsda8QFdW2NaNeunTp06CDJsWuE9aBq7969ioiIKFFv/X/+/fffq2nTpqwRqltrhFVNH0ccP35cUsWh19qONeKGMuv/im1gAIA6KyMjQ9nZ2WVu4Zk8ebJOnjypjIwMJ40QRhcaGqqAgAC9+OKLJT69tlgsevHFF3X55ZcrMDDQSSNEXdCsWTNJKvMDDWu5tR0qz1RU3jkz2Pz8889lnrbi9O0FjRs31tmzZx3er2SMbRts8ah628pu8fjvf/+rtm3bqrCwUPv27dNll13GNrBqtDX6GlETbb/44guNHz9eubm58vb2LlGfm5urRo0aac6cObZPmFkjjNPWCNvAXFxclJycrLi4OEVHR2vSpEm2rYSzZs1SQkKCbSsha0TV20qXxjawwsJC9ejRQ+3atSt2lle6MLdjYmL0ww8/aPPmzbZtz5f6GmGxWOTl5aUrrriizLYS28Aqzd3dvdJ77C7+43Fk24sXBiO2/etrqakxuLq62t4s6ltbFxeXSu/jrE9tCwoKii1sJpOp0n8bRmgr1dzffX1aI2qibfPmzSVVvIWnefPmpf4ujfB3zxpR9baO/lvu06ePFi1apBkzZti2jUpSy5YttWjRIvXp08euMUisEUZqWxN/n66urpo6dari4uIUExNTauhdtGiRrS/WiMrjzEol/frrr1W+IOhS07hx4zL3twLlyc/PV+vWrSVd2G/r6enp5BGhrqjqp5lAeXJzc3XjjTdKkj744AP16NGDeYMqSUpK0owZM4pda9OyZUtNmTKl1NB7KavsmRWOvAEAdZb108yEhATFxMQUe7ZBTEyMEhISNGXKFA44USkXz5MuXbowb1Blffr00YYNG2xff/DBB9q8eTNBxQ5sAwMA1GnV2cIDADWF0OtYhBUAQJ3Xp08fRUZGsoUHAOoZtoEBAOoFPs0EgPqHsAIAAADAkAgrAJzu4tsWb9++vdz7swMAgEsHYQWAUyUlJal37962r0eOHKkePXooKSnJiaMCAABGQFgB4DRJSUmKi4tT+/bti91ytl27doqLiyOwAABwiSOsAHCKwsJCzZw5U9HR0VqzZo0iIiLk7e2tiIgIrVmzRtHR0XrhhRfYEgYAwCWMsALAKdLS0nTo0CFNnjy52FPHJcnFxUWTJk1SZmam0tLSnDRCAADgbIQVAE5x/PhxSVJYWFip9dZyazsAAHDpIawAcIpmzZpJkvbu3VtqvbXc2g4AAFx6eIJ9FVgsFmcPwdAKCwv5HaHSOnXqpICAAL344otas2ZNsa1gFotFs2bNUsuWLdWpUyfmFSrl4nlisViYN6gy5hAcgXlUOZX9vRBWKsnFxUUWi4WLfcthNpv5/aBKJk2apLFjxyomJkaTJk1SWFiY9u7dq1mzZikhIUELFy6UJOYVKuXieVJYWMi8QZU1aNBABw4cUKNGjZSTk8McQrWwFlWeq6trhW0IK5Xk4+OjoqIiZw/D0Pz8/Jw9BNQxgwYN0mWXXaZnnnlG3bp1s5UHBgbqvffeU3R0tBNHh7rm4jc9Hx8feXl5OXE0qMt4P4M9WIsqz9fXt8I2piKOwCvl1KlTzh6C4fn5+fF7QrXk5OQoKChIkrRy5Ur16tWrUp+2ABc7e/asAgICJEmHDh3iAAHVxvsZHIF5VLHKfDDABfYAnO7iYNK1a1eCCgAAkERYAQAAAGBQhBXY7ezZs2rSpIlMJpPOnj3r7OEAAACgniCsAAAAADAkwgoAAAAAQyKsAAAAADAkwgoAoF7w8vJSdna2ioqKuG0xANQTdeKhkLt379b8+fO1a9cumc1mBQcHa9SoUYqKiqp0H8eOHdM777yjb7/9VkeOHJGnp6cCAwN19913q1+/ftwqFQAAADAYw4eV1NRUxcbGysPDQ3379pWXl5eSk5M1fvx4HT16VGPGjKmwj0OHDmnQoEE6ffq0IiMj1atXL+Xl5enf//63Jk6cqO3bt2vWrFm18GoAAAAAVJahn2BvNpvVp08fHT16VCtXrlRoaKgkKTc3VwMHDlRWVpY2bNigFi1alNvPc889p08++USTJ0/WyJEjbeU5OTnq37+/jhw5oo0bN5bbD08gLRtPjYa9mENwJJ4aDXsxh+AIzKOK1fkn2KempiozM1PR0dG2oCJJPj4+iouLU0FBgeLj4yvs59ChQ5Kknj17Fitv1KiROnToIIkwAgAAABiNocNKWlqaJCkyMrJEnbUsPT29wn6Cg4MlSZs3by5WnpOTo127dqlp06a6/vrr7R0uAAAAAAcy9DUrBw8elCQFBgaWqGvatKk8PT2VkZFRYT/33nuvNm7cqFmzZmnLli0KCQmxXbPSsGFDLViwQA0bNnT08AEAAADYwdBhJS8vT9KFbV+l8fb2Vm5uboX9XHHFFfr000/15JNP6ptvvtGWLVskSQ0bNtSQIUPUunXrCvvw9fWVi4uhT0Q5jYeHh+3ffn5+XG+AKmMOwdEqsw8aKA9zCI7APLKfocOKo2RkZCguLk6enp5atmyZQkNDlZubqy+++EJvvPGGUlJStGzZsnJvX3zmzJlaHHHdcvbsWdu/T506pfPnzztxNKiLmENwJC5qhb2YQ3AE5lHFKhPmDB1WvL29JanMsyd5eXny9fWtsJ+nn35aR44c0VdffaWmTZtKuvDwsPvvv19//PGHPvjgA61fv1533HGH4wYPAAAAwC6G3tcUFBQkSaVel3LixAnl5+eXej3LxfLy8rRz505dd911tqBysS5dukiSfvrpJ/sHDAAAAMBhDB1WwsPDJUkpKSkl6qxl1jZlKSgokFT2rYmzs7MlFd8zDwAAAMD5DB1WunbtqoCAACUkJBQ785Gbm6tFixbJ3d1dMTExtvLjx4/rwIEDxbaN+fn56ZprrtGRI0e0atWqYv3n5OTovffek/S/MywAAAAAjMHQ16y4ublp5syZio2N1bBhw9S3b195eXkpOTlZWVlZmjhxovz9/W3tZ8+erfj4eM2aNUsDBgywlU+aNEkPPfSQpk6dqvXr1ys0NFQ5OTnauHGjsrOz1bt3b3Xr1s0ZLxEAAABAGQwdViQpIiJCy5cv17x585SYmCiz2azg4GA98cQTioqKqlQfPXv21CeffKIlS5Zox44dSk9Pl4eHh6677jqNHTtW99xzTw2/CgAAAABVZSoqKipy9iDqAm49V7azZ88qICBAknTo0CGekYEqYw7BkbhdKOzFHIIjMI8qVplbFxv6mhUAlwYvLy9lZ2erqKiIoAIAAGwIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACuxUWFtr+vW3btmJfAwAAANVFWIFd1q1bp+7du9u+Hjx4sDp27Kh169Y5cVQAAACoDwgrqLZ169Zp1KhRateunbZt26bc3Fxt27ZNbdu21ahRowgsAAAAsAtPsK8knkBaXGFhoTp27Ki2bdtqzZo1cnH5X+61WCyKiYnR7t279d1338nV1dWJI0VdwtN+4QjMI9iLOQRHYB5VjCfYo8Zs27ZNmZmZmjx5crGgIkkuLi6aNGmSMjIytG3bNieNEAAAAHUdYQXVcuzYMUlSWFhYqfXWcms7AAAAoKoIK6iWK6+8UpK0d+/eUuut5dZ2AAAAQFURVlAtXbt2VcuWLfXiiy/KYrEUq7NYLJo1a5YCAwPVtWtXJ40QAAAAdR1hBdXi6uqqGTNmKCEhQTExMcXuBhYTE6OEhARNnz6di+sBAABQbW7OHgDqrn79+mnp0qWaOnWqunXrZisPDAzU0qVL1a9fPyeODgAAAHUdYQV26devn3r27KmgoCBJ0sqVK9WrVy/OqAAAAMBubAOD3S4OJl27diWoAAAAwCEIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMibACAAAAwJAIKwAAAAAMyc3ZA6iM3bt3a/78+dq1a5fMZrOCg4M1atQoRUVFVamfkydP6u2339bXX3+t33//XZ6engoKClL//v01dOjQGho9AAAAgOowfFhJTU1VbGysPDw81LdvX3l5eSk5OVnjx4/X0aNHNWbMmEr189NPP2nMmDHKyclRz5491bt3b+Xn5+vAgQPatGkTYQUAAAAwGEOHFbPZrGnTpslkMmnZsmUKDQ2VJI0dO1YDBw7U7Nmz1bt3b7Vo0aLcfvLy8vTQQw9JklavXq3WrVuX+DkAAAAAjMXQ16ykpqYqMzNT0dHRtqAiST4+PoqLi1NBQYHi4+Mr7Gf58uU6cuSIHn/88RJBRZLc3Ayd2QAAAIBLkqGP0tPS0iRJkZGRJeqsZenp6RX2k5iYKJPJpN69e+vXX3/V1q1bde7cOV177bXq3r27PDw8HDtwAAAAAHYzdFg5ePCgJCkwMLBEXdOmTeXp6amMjIxy+zh//rx++eUXNWnSRB999JHmz58vi8Viqw8ICNDChQsVEhLi0LEDAAAAsI+hw0peXp6kC9u+SuPt7a3c3Nxy+zhz5owKCwt1+vRpvfnmm3ryySfVv39/mc1mrVixQm+99ZYefPBBJSUlqUGDBmX24+vrKxcXQ++ac5qLz0z5+fnJy8vLiaNBXefn5+fsIaAeYB7BXswhOALzyH6GDiuOYD2LUlhYqGHDhhW7e9ijjz6q3377TUlJSfrXv/6l/v37l9nPmTNnanysddXZs2dt/z516pTOnz/vxNGgLvPz89OpU6ecPQzUccwj2Is5BEdgHlWsMmHO0KcKvL29JanMsyd5eXllnnWxurj+lltuKVFvLdu7d291hwkAAACgBhg6rAQFBUlSqdelnDhxQvn5+aVez3IxT09PXXnllZKkRo0alai3lv355592jhYAAACAIxk6rISHh0uSUlJSStRZy6xtyhMRESFJ2r9/f4k6a1lFz2oBAAAAULsMHVa6du2qgIAAJSQk6KeffrKV5+bmatGiRXJ3d1dMTIyt/Pjx4zpw4ECJbWNDhgyRJL3zzjvKycmxlZ84cUIffvihXFxcdPvtt9fsiwEAAABQJYa+wN7NzU0zZ85UbGyshg0bpr59+8rLy0vJycnKysrSxIkT5e/vb2s/e/ZsxcfHa9asWRowYICtvEOHDho9erTef/993XHHHerVq5fMZrP+/e9/6+TJk5owYYKuueYaZ7zEesHLy0vZ2dlcSAYAAACHqnJY+evzSEwmk7y8vHTdddcpKipKw4YNk7u7e7E2iYmJWrp0qQ4fPiyz2axWrVrp8ccfV4cOHSr8eREREVq+fLnmzZunxMREmc1mBQcH64knnlBUVFSlx/30008rODhYy5YtU3x8vEwmk0JDQ/X888/rtttuq3Q/AAAAAGqHqaioqKgq32ANK3feeaekC7cEzsrK0q5du2SxWNS1a1e9++67cnP7Xw5KTEyUn5+fOnXqpMLCQj3//PPasGGDkpKSbBe/Gx1nDCrGmRXYizkER2AewV7MITgC86hilbl1cbW3gb300kvFvv7hhx80fPhwbdu2TevXry/2zJKLz4C4u7tr5MiR+vzzz7Vnz546E1YAAAAA1C6HXWDftm1b29mW0u7edbENGzbIzc2txJYyAAAAALBy6N3AWrVqJUnKzs4us813332nd955RyNGjFBAQIAjfzwAAACAesShYeXs2bOSpCZNmpRaf+DAAY0dO1bh4eF6/PHHHfmjAQAAANQzDg0rW7ZskSR17969RF1mZqZGjRqlVq1aaeHChcUuwAcAAACAv7I7MVgsFh0+fFhLlixRenq6/v73v5e4pfChQ4c0YsQIXX/99Vq4cKE8PT3t/bEAAAAA6rlqh5XSLo4fPHiwpk+fLpPJVKz8zTff1O+//66zZ8+qV69etvLY2Fjdd9991R0CAAAAgHrM7ues/Pnnn9q3b59+/fVXSSrx9Pj6gvtkV4z7icNezCE4AvMI9mIOwRGYRxWr1eesvPvuu3r11Vc1ffp0denSRS1atKhu1wAAAADguAvsY2NjFRkZqf/+979asGCBo7oFAAAAcIly6N3ArLcj/uKLL5SVleXIrgEAAABcYhwaVm644QbdeuutMpvNevfddx3ZNQAAAIBLjEPDiiQ9/PDDMplMWr16tU6cOOHo7gEAAABcIhweVlq3bq3bbrtNf/75p95//31Hdw8AAADgElHlu4H9/PPPFbaZP39+tQYDAAAAAFYOP7MCAAAAAI5AWAEAAABgSIQVAAAAAIZEWAEAAABgSIQVAAAAAIZEWAEAAABgSIQVAAAAAIZEWAEAAABgSIQVAAAAAIZEWAEAAABgSIQVAAAAAIZEWAEAAABgSIQVAAAAAIZEWAEAAABgSIQVAAAAAIZEWAEAAABgSIQVAAAAAIZEWAEAAABgSIQVAAAAAIZEWAEAAABgSIQVAAAAAIZEWAEAAABgSHUirOzevVv33XefOnXqpHbt2mnw4MFKTEysdn9nzpxR9+7dFRISonvvvdeBIwUAAADgKG7OHkBFUlNTFRsbKw8PD/Xt21deXl5KTk7W+PHjdfToUY0ZM6bKfU6fPl15eXk1MFoAAAAAjmLoMytms1nTpk2TyWTSsmXLNGPGDD399NNau3atgoKCNHv2bGVlZVWpzw0bNighIUFPPPFEDY0aAAAAgCMYOqykpqYqMzNT0dHRCg0NtZX7+PgoLi5OBQUFio+Pr3R/2dnZeu6559S/f3/17NmzJoYMAAAAwEEMHVbS0tIkSZGRkSXqrGXp6emV7u/ZZ5+Vq6urpkyZ4pgBAgAAAKgxhr5m5eDBg5KkwMDAEnVNmzaVp6enMjIyKtXX2rVrlZycrIULF8rX11e5ubmOHCoAAAAABzN0WLFeBO/j41Nqvbe3d6VCx7Fjx/TCCy8oOjpat956a7XG4uvrKxcXQ5+IMgQ/Pz9nDwF1HHMIjsA8gr2YQ3AE5pH9DB1WHGXq1Klyc3Oza/vXmTNnHDii+snPz0+nTp1y9jBQhzGH4AjMI9iLOQRHYB5VrDJhztBhxdvbW5LKPHuSl5cnX1/fcvuIj4/XN998o7lz56pJkyYOHyMAAACAmmHosBIUFCRJysjIUFhYWLG6EydOKD8/X23atCm3jx9//FGS9Oijj5Zan5KSopCQELVu3Vpr1661f9AAAAAAHMLQYSU8PFxvv/22UlJS1Ldv32J1KSkptjblad++vfLz80uU5+fnKzExUVdddZUiIyN19dVXO27gAAAAAOxm6LDStWtXBQQEKCEhQSNGjLA9ayU3N1eLFi2Su7u7YmJibO2PHz+u3NxcNWvWzHZRflRUlKKiokr0ffjwYSUmJur666/XCy+8UCuvBwAAAEDlGfr2Vm5ubpo5c6aKioo0bNgwTZs2TS+99JL69++vgwcPasKECfL397e1nz17tqKiovTll186cdQAAAAAHMHQZ1YkKSIiQsuXL9e8efOUmJgos9ms4OBgPfHEE6WeMQEAAABQP5iKioqKnD2IuoBbz1WMW/TBXswhOALzCPZiDsERmEcVq8ytiw29DQwAAADApYuwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADImwAgAAAMCQCCsAAAAADMnN2QOojN27d2v+/PnatWuXzGazgoODNWrUKEVFRVX4vUVFRfrmm2+0ceNG7dy5U0eOHJHZbFZgYKCioqI0evRoNWjQoBZeBQAAAICqMHxYSU1NVWxsrDw8PNS3b195eXkpOTlZ48eP19GjRzVmzJhyv//8+fO6//775eHhoc6dOysyMlLnz59XSkqK5syZo6+++kofffSRLrvsslp6RQAAAAAqw9BhxWw2a9q0aTKZTFq2bJlCQ0MlSWPHjtXAgQM1e/Zs9e7dWy1atCizDxcXFz322GMaOnSofH19beUFBQV6+OGHtWnTJi1btkyxsbE1/noAAAAAVJ6hr1lJTU1VZmamoqOjbUFFknx8fBQXF6eCggLFx8eX24e7u7sefPDBYkHFWv7AAw9IktLT0x0/eAAAAAB2MXRYSUtLkyRFRkaWqLOW2RM03NwunFhydXWtdh8AAAAAaoahw8rBgwclSYGBgSXqmjZtKk9PT2VkZFS7/9WrV0uSbr755mr3AQAAAKBmGPqalby8PEkXtn2VxtvbW7m5udXqe/Pmzfr000913XXXadCgQRW29/X1lYuLobOdIfj5+Tl7CKjjmENwBOYR7MUcgiMwj+xn6LBSU3bv3q3x48fLx8dHc+fOlYeHR4Xfc+bMmVoYWd3m5+enU6dOOXsYqMOYQ3AE5hHsxRyCIzCPKlaZMGfoUwXe3t6SVObZk7y8vDLPupRlz549uvfee+Xi4qJ3331XrVq1snucAAAAABzP0GElKChIkkq9LuXEiRPKz88v9XqWsuzZs0djxoyRxWLRkiVL1KZNG0cNFQAAAICDGTqshIeHS5JSUlJK1FnLrG0qYg0qhYWFevfdd9W2bVvHDRQAAACAwxk6rHTt2lUBAQFKSEjQTz/9ZCvPzc3VokWL5O7urpiYGFv58ePHdeDAgRLbxvbu3asxY8bIbDbrnXfeUfv27WvrJQAAAACoJkNfYO/m5qaZM2cqNjZWw4YNU9++feXl5aXk5GRlZWVp4sSJ8vf3t7WfPXu24uPjNWvWLA0YMECSdPr0aY0ZM0Y5OTnq3r27vv32W3377bfFfo6Pj49GjRpVmy8NAAAAQAUMHVYkKSIiQsuXL9e8efOUmJgos9ms4OBgPfHEE4qKiqrw+/Py8mx38tqyZYu2bNlSok2LFi0IKwAAAIDBmIqKioqcPYi6gFvPVYxb9MFezCE4AvMI9mIOwRGYRxWr87cuBgAAAHDpIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDIqwAAAAAMCTCCgAAAABDcnP2ACpj9+7dmj9/vnbt2iWz2azg4GCNGjVKUVFRle7j/PnzWrx4sb744gv9/vvv8vX1Va9evfTYY4/p8ssvr8HRAwAAAKgOw4eV1NRUxcbGysPDQ3379pWXl5eSk5M1fvx4HT16VGPGjKmwD4vFogcffFApKSlq166dbr/9dmVkZGjVqlXatm2bVq5cqSZNmtTCqwEAAABQWYYOK2azWdOmTZPJZNKyZcsUGhoqSRo7dqwGDhyo2bNnq3fv3mrRokW5/cTHxyslJUXR0dF67bXXZDKZJEmffPKJnnvuOb3xxhuaPn16jb8eAAAAAJVn6GtWUlNTlZmZqejoaFtQkSQfHx/FxcWpoKBA8fHxFfazatUqSdKECRNsQUWShgwZooCAAK1bt07nzp1z/AsAAAAAUG2GDitpaWmSpMjIyBJ11rL09PRy+/jzzz/1ww8/6JprrilxBsZkMqlbt27Kz8/X3r17HTRqAAAAAI5g6LBy8OBBSVJgYGCJuqZNm8rT01MZGRnl9pGZmSmLxaKgoKBS663l1p8FAAAAwBgMfc1KXl6epAvbvkrj7e2t3Nzccvuw1nt7e5fZx8U/qyx+fn7l1uMCfk+wF3MIjsA8gr2YQ3AE5pH9DH1mBQAAAMCly9BhxXrWo6yzJ3l5eWWedbGy1pd15sRaXtaZFwAAAADOYeiwYr2epLTrUk6cOKH8/PxSr2e5WEBAgFxcXMq8JsVaXtY1LQAAAACcw9BhJTw8XJKUkpJSos5aZm1TloYNG6pNmzb67bfflJWVVayuqKhI3377rTw9PRUWFuagUQMAAABwBEOHla5duyogIEAJCQn66aefbOW5ublatGiR3N3dFRMTYys/fvy4Dhw4UGLb2ODBgyVJs2fPVlFRka18xYoVOnTokPr166eGDRvW7IsBAAAAUCWGDitubm6aOXOmioqKNGzYME2bNk0vvfSS+vfvr4MHD2rChAny9/e3tZ89e7aioqL05ZdfFuvnzjvvVGRkpBISEjRkyBC99tpreuSRR/T888/L399fjz32WC2/svpj+PDhCgkJKVa2fft2hYSEaP78+U4aFeqj+fPnKyQkRNu3b3f2UFAHMF9Q0xz1XsdcRUUu9WMtQ4cVSYqIiNDy5cvVoUMHJSYm6pNPPtHll1+uOXPmaMyYMZXqw8XFRW+99ZYefvhhZWdna+nSpdq5c6cGDhyoTz/9VE2aNKnhV1H7Dh8+rJCQEN17771ltrFO9GeeeaYWR4ZLxeLFixUSEqLFixeXWt+/f/9y599DDz2kkJAQbdu2rSaHCScpKipScnKyxo0bpx49eigsLEzt27fXHXfcoRdffFH79+939hBRB1yK68yldJBal1j/v1T2v+HDhzt7yHWGoZ+zYtWmTRu9++67FbZ76aWX9NJLL5Va5+HhoXHjxmncuHGOHh6AUnTp0kXShQX8/vvvL1Z36tQp/fzzzzKZTEpLSyvxvRaLRd999508PDzUoUMHSdKwYcMUFRWl5s2b1/zgUaNOnz6tRx99VKmpqWrUqJG6deumgIAAFRQUaP/+/Vq+fLk++ugjLV261DaPgNI4ep2pijZt2igxMZHnaECS1KJFixLHmDk5Ofrwww/VokUL3XnnnSXaV9bLL7+s//73vw4ZZ11UJ8IKgNr39NNPKz4+Xj///HO1vv/GG2+Ul5eXdu7cKbPZLDe3/y036enpKioq0u23367k5GQdP35czZo1s9Xv27dPZ86cUefOndWgQQNJUpMmTerlWdC6xN45IUlms1njxo1Tenq67rjjDj377LMlbh1//PhxzZkzp8KH/qLuM9o6UxWXXXaZrrvuumqNG8Zkz3z09/fXww8/XKzs8OHDtrDy17qquNQ/pDP8NjDUrr1792r69OmKjo5Wx44d1aZNG/Xr10+LFy9WQUGBXX2Xd9rzlltu0S233GL7OiMjQ+3bt1ePHj106tSpYm3Lq4NxuLm5qVOnTsrPz9eePXuK1aWlpalhw4a67777JKnEXm3rp6AXf6pe2r5u63bHp59+WhkZGRo7dqzCw8PVrl07jRo1Svv27St1bGlpaRo2bJjatWunLl266LHHHtPvv/9e6r5gONbatWuVnp6u8PBwvfzyy6U+46pZs2aaNWuWevToYSuryvpRlovny4EDB/TAAw+oU6dOCg8P14QJE5SdnS1J2rVrl0aOHKkOHTooPDxcU6ZMUX5+fql9pqenKy4uTl26dFFYWJhuv/12zZkzp8SnoBdv3dmzZ49Gjx6t9u3bq2PHjho7dqwOHz5c4fhRkqPXGUn66quvNHLkSIWHh+umm25SdHS0lixZosLCwmLtytuOVd01Zt26derfv7/atGmjyMhIzZw5U+fOnbPVz58/XyNGjJAkLViwoNi2IuZQ3VHVY61L/b2JsIJiVq5cqS+//FLBwcG6++67NXDgQBUVFen111/XhAkTam0cgYGBmjZtmo4dO6apU6faygsKCjRhwgSdO3dOr7zyCqffDe7iLRoX2759u9q2bauwsDD5+vqWqE9NTS32/RXJysrS4MGDdebMGd11113q1q2btm3bphEjRuiPP/4o1jYlJUWjR4/W7t271bt3bw0ePFhHjhzR0KFDlZOTU92Xikr67LPPJEkPPvigXFzKfwvy8PCokTEcPnxYQ4YM0fnz5zVo0CC1bt1a69ev19ixY/Xdd99p1KhR8vT01N13362AgAB99tlnmjFjRol+li9fruHDh2vnzp3629/+puHDh+vKK6/UokWLNHr0aJ0/f77E9+zZs0f//Oc/5e7uriFDhigsLExfffWVRo8erT///LNGXm9958h15vXXX9fYsWP122+/6bbbbtPQoUPVoEEDvfLKKxo/fnylxlPdNWbZsmV65plndP311+uee+5Ro0aN9NFHH2nKlCm2Np07d7ZtJ+rcubNte/u4cePUqFGjSo0PzmeUY626gm1g9VxmZmaZF+H99bkzkhQXF6dnn31Wrq6utrKioiJNmTJFq1ev1o4dO9SxY8caG+/FBgwYoJSUFK1fv17Lly/X0KFDNWfOHO3du1cPPPCAIiIiamUcqD7rQUBaWpri4uIkSdnZ2frPf/6jcePGycXFRR07dix2EGGxWLRjxw41bNhQbdu2rdTPSUtL0+OPP15sz/obb7yht956S59//rmtvLCwUM8884wKCwv1wQcfqFOnTrb2EydO1Jo1a+x9ySiH2WzWnj17bJ+GO0t6eromT56skSNHSrqwxj3wwAPavHmzHnzwQb3++uu69dZbJV34gOSuu+7SF198occff1xXXHGFJGn//v164YUXFBISoqVLlxb74GTx4sV6/fXX9fHHH5e4EczmzZs1Z84cRUVF2cqeeuoprV27Vl999ZX69u1b0y+/3nHUOrN161YtXrxYkZGRmj9/vjw9PSVdmB/PPfecVqxYoQ0bNqh3795ljsWeNebbb7/V6tWrde2110qSxo8fr/79+ysxMVFPPfWUrrzySttrjY+PV+fOne3aWgTnMdKxVl3AmZV6LjMzUwsWLCj1v/j4+BLtmzdvXuyPR5JMJpOGDRsmSbV+x5Tnn39eLVq00Msvv6yPPvpI7733ntq0aaNHHnmkVseB6rnhhhvUqFEj7dy50/Ypc1pamoqKimxvup07d1ZmZqZ+//13SdJPP/2knJwctW/fvtKfrPv7+ys2NrZY2cCBAyWp2NaQHTt2KCsrS7169SpxsPzYY4+VmPtwrNOnT6ugoEB+fn7VukbAUVq2bGnbSiNdWOOs4SE0NNQWVCTJ3d1dvXv3ltlsLnaHshUrVshsNmvatGklzvDGxsaqSZMmSkhIKPGzw8PDiwUVSbrrrrskqcQ2JlSOo9aZjz/+WJI0Y8YMW1CRLsyPJ554QiaTSevXry93LPasMSNGjLAFFenCQ62jo6NlsVj0f//3f5X9daAOMNqxltFxZqWei4yM1JIlS0qt2759e7E3bEk6f/68li1bpvXr1+vXX39Vfn5+sQdpHj9+vEbH+1c+Pj567bXX9M9//lMzZ86Ul5eXXn/99WIXUcJ+t9xyS6ln2iSVuk921qxZGjBgQIX9Wj/R3LRpk3bv3q1OnTpp+/btatCgge3TzPDwcEkX5mNMTEyZ+8jLExoaWmJL0VVXXSVJxbZdWK9hKe0Tq6uvvlpXX301+77/v5qaE0YQEhIik8lUrMx64XVoaGiJ9ta6i9e/H374QZK0ZcuWUg8s3Nzc9Ntvv5Uov/HGG0uUlTZX6yOjrzM//PCDPD09tXr16lJ/TsOGDfXrr7+WOxZ71phLeW44gzPXOKMdaxkdR3wo5pFHHtGmTZsUFBSkqKgoXX755XJzc7Pdfq+0Pdg17cYbb1Tz5s116NAh9ejRQy1btqz1MdR3I0aMKHHnpa+++kr79u0r9XbfpR3QlSUiIkKbNm1SWlqaOnXqpLS0NLVt29b2aWZoaKh8fHxsBxHWrRpV2eZX2kXa1kBrsVhsZXl5eZKkyy+/vNR+rrjiCsLK/1cTc6Jx48Zyd3fX6dOndf78+Rq7JqUipc0X66ec5dWZzWZb2ZkzZyRJixYtctjPvniu1kdGX2fOnDkjs9msBQsWlPlzyrrRgpU9a8ylPDecoSbnY0WMeKxlZIQV2OzevVubNm1SZGSkFi9eXOwU5ffff68PP/zQrv5NJlOxN/uL5ebmysfHp9S6V155RYcOHVLjxo2VlJSkO++8Uz179rRrLChu1KhRJcqysrK0b98+u/dEX3zx69133639+/cX69PV1dW2n9z63ANPT0/ddNNNdv3c0lgPBk6ePFlq/V8vxr+U1cSccHNz00033aSdO3cqPT1dN998c6W/t7rrR02xzqUdO3aUepCJkoy+zlj/P9rzJHnWmLqjJudjeWr6WKs+4poV2Bw6dEiS9Le//a3EXsrvvvvO7v59fX117NixEuWHDx8u8xT3119/rY8//lidO3fW6tWr5evrq0mTJrHg1yEhISFq3Lixvv/+e23ZskXShf3jFwsPD1dWVpaSk5OVm5urjh071shWv9atW0uSdu7cWaLu6NGjtv3sqDnWa4kWLVpUbNtDaS7+dLE660dNatOmjaT/bQeDczlinWnTpo1Onz6tgwcPVnsctbHGWN+f/3orZdQNNX2sVR8RVmBjfejQjh07ipX/5z//0eLFi+3uPywsTFlZWcWeJHz+/Hm99NJLpbY/ceKEJk2aJF9fX7366qvy9/fX9OnTdfLkSU2cOLHCAx0Yg4uLi8LDw3Xu3DktWbJEDRo0ULt27Yq1sR5ULFy4UFLVrlepio4dO6p58+batGmTdu3aVaxu7ty5vPnXgv79+9u26UyaNMm2beZif/zxh6ZOnapvvvnGVlbV9aOmDR06VG5ubpoxY4aOHDlSoj4nJ0c//vijE0Z2aXLEOmN9js/kyZNLfYbXiRMndODAgXLHURtrjK+vr6QL4Qd1T00fa9VHbAODTZs2bdSmTRslJSXpxIkTatu2rX7//Xdt3LhRPXv21IYNG+zqf/To0dq6davuv/9+9e3bV5dddpm2bt2qRo0aqWnTpsXaFhUVaeLEicrOzta8efNsFxn+4x//0MCBA/XZZ5/p/fffL3FbUBhTly5d9OWXX+qXX35R586dS1yrcMMNN8jT01O//PKLrX1NcHV11XPPPaeHHnpII0eOVFRUlJo2bar09HQdO3ZMrVu3tuvp7KiYm5ubFi5cqEcffVTx8fHauHGjbr75Zvn7+6ugoED79+9XWlqazGaz7rjjDtv3VWX9qA3BwcF69tln9dxzz+kf//iHevbsqYCAAJ09e1aHDx9WWlqa7rzzTk2fPr3Wx3apsned6dGjhx566CG9+eabuv3229W9e3c1b95cp0+fVkZGhnbs2KHHHnus3KfW18Yac+2116pZs2Zav369PDw8dOWVV8pkMmn48OG1vh0SVVfTx1r1EWdWYOPq6qq3335bd911lzIzM/Xxxx9r//79euqpp/Tkk0/a3X9kZKTeeOMNBQQEaO3atfrXv/6lm2++We+9916JN5X33ntPW7du1aBBg0rc037KlCkKCgrS7Nmz+eSyjrj4oKC0IOLm5qYOHTpIurDnu7S74jhKz549tWTJEoWFhSkpKUkrV67UlVdeqeXLl8tisXD9QS1o3Lixli5dqnnz5ik8PFzfffed3n//fa1YsULHjh3T4MGDtW7dumLbeKqyftSWwYMHa8WKFbr11ltte803bNigU6dOadSoUbbnuKB2OGKdefTRR/X++++rY8eO2rZtm5YuXaqvv/5aBQUFGjdunPr161fhOGp6jXF1ddWCBQvUrl07JSQkaN68eZo7d67tpg8wtpo+1qqPTEXspQEA5eXl6eabb1ZwcLBWrVrl7OEAqGdYY4Dq4cwKgEtKfn5+ieskCgsL9corr+jcuXPFHggIAFXFGgM4FtesALikZGRkaOjQoYqMjJS/v7/Onj2rHTt2aP/+/WrVqpXtIlsAqA7WGMCx2AYG4JKSnZ2tV199VWlpaTp58qTMZrOaN2+uv//973rwwQfVqFEjZw8RQB3GGgM4FmEFAAAAgCFxzQoAAAAAQyKsAAAAADAkwgoAAAAAQyKsAAAAADAkwgoAAAAAQyKsAAAAADAkwgoAAAAAQyKsAAAAADCk/weJo6wHLfNWOAAAAABJRU5ErkJggg==", "text/plain": [ "
<xarray.DataArray 'y' ()> Size: 8B\n", "array(97.39806397)
<xarray.DataArray 'y' ()> Size: 8B\n", "array(97.39806397)