Hierarchical Regression With Missing Data

hierarchical-regression-with-missing-data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import numpyro\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import jax.numpy as jnp\n",
    "import numpyro.distributions as dist\n",
    "\n",
    "from jax import random\n",
    "from numpyro.infer import Predictive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 500\n",
    "\n",
    "# Generate machine process parameters for both groups\n",
    "feed_speed = np.random.uniform(1, 10, n_samples)\n",
    "pull_acceleration = np.random.uniform(0.5, 5, n_samples)\n",
    "cutting_wait_time = np.random.uniform(0.1, 2, n_samples)\n",
    "\n",
    "# Create product groups\n",
    "product_group = np.array(['Product_Group_1'] * (n_samples // 2) + ['Product_Group_2'] * (n_samples // 2))\n",
    "\n",
    "# Calculate production speed\n",
    "production_speed = np.zeros(n_samples)\n",
    "\n",
    "# Product Group 1: depends on feed speed and pull acceleration (Cutting Wait Time is NaN)\n",
    "mask_pg1 = product_group == 'Product_Group_1'\n",
    "\n",
    "# Increased coefficients for Product Group 1 to ensure higher production speed\n",
    "production_speed[mask_pg1] = (feed_speed[mask_pg1] * 2.5 + \n",
    "                              pull_acceleration[mask_pg1] * 3 + \n",
    "                              np.random.normal(0, 0.5, sum(mask_pg1)))\n",
    "\n",
    "# Product Group 2: depends on all three parameters\n",
    "mask_pg2 = product_group == 'Product_Group_2'\n",
    "production_speed[mask_pg2] = (feed_speed[mask_pg2] * 1 + \n",
    "                              pull_acceleration[mask_pg2] * 0.5 + \n",
    "                              cutting_wait_time[mask_pg2] * 3 + \n",
    "                              np.random.normal(0, 0.5, sum(mask_pg2)))\n",
    "\n",
    "# Create DataFrame with NaN for Cutting Wait Time in Product Group 1\n",
    "simulated_dataset = pd.DataFrame({\n",
    "    'Product_Group': product_group,\n",
    "    'Feed_Speed': feed_speed,\n",
    "    'Pull_Acceleration': pull_acceleration,\n",
    "    'Cutting_Wait_Time': [np.nan if x == 'Product_Group_1' else cutting_wait_time[i] for i,x in enumerate(product_group)],\n",
    "    'Production_Speed': production_speed\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAE8CAYAAACIDWb/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMDklEQVR4nO3deVhUZf8/8PeAMOzDKjCyiKAoGtijZmhuSCKVaVGplYH5mBbiQmpRmVvmVkaL6fdpkTZb7HHLTFQUzDXFUHEh4UExAVFkRwaE8/vDHydHQJlhZs4A79d1zRVnuc/5nHuO8Omce5EJgiCAiIiIiAzKROoAiIiIiNojJmFEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEEuvcuTOioqIkOfeCBQsgk8kkOXdTLly4AJlMhoSEBKlDMQqG/I6GDh2KoUOHisvJycmQyWT4+eefDXL+qKgodO7c2SDnIjIGTMKI9OTUqVN46qmn4O3tDQsLC3Tq1AkPP/wwPv74Y6lD04lPP/1U0kSpPkGo/8jlcri6umLo0KF49913cfXqVa2PfebMGSxYsAAXLlzQXcAAEhIS1GK2sLCAUqlEWFgYPvroI5SVlenkPLm5uViwYAHS0tJ0cjxdMubYiAytg9QBELVFBw8exLBhw+Dl5YXJkyfDzc0Nly5dwuHDh/Hhhx8iJiZG3DcjIwMmJq3v/4c+/fRTODs7S/YUr9706dPRr18/1NbW4urVqzh48CDmz5+PVatW4aeffkJISIjGxzxz5gwWLlyIoUOH6uXJzKJFi+Dj44Oamhrk5+cjOTkZM2fOxKpVq7B161YEBgaK+7711lt4/fXXNTp+bm4uFi5ciM6dO6N3797NLrdz506NzqONu8X22Wefoa6uTu8xEBkLJmFEerBkyRIoFAocPXoU9vb2atsKCgrUluVyuQEja3sGDRqEp556Sm3diRMnMGLECERERODMmTNwd3eXKLrGhYeHo2/fvuJyXFwc9uzZg8ceewyPP/44zp49C0tLSwBAhw4d0KGDfn9VV1ZWwsrKCubm5no9z72YmZlJen4iQ2t9//tN1ApkZWWhZ8+eDRIwAOjYsaPa8p1twupfWe3fvx/Tp0+Hi4sL7O3tMWXKFFRXV6O4uBgvvPACHBwc4ODggLlz50IQBLF8/Wu65ORktfM0t63VunXrEBISgo4dO0IulyMgIABr1qxpEPPp06eRkpIivlq7vS1RcXExZs6cCU9PT8jlcvj5+WH58uUNnnIUFxcjKioKCoUC9vb2iIyMRHFx8V3ja46goCDEx8ejuLgYn3zyibj+4sWLeOWVV+Dv7w9LS0s4OTnh6aefVnvtmJCQgKeffhoAMGzYMPH66utzy5YtePTRR6FUKiGXy+Hr64vFixejtra2RTGHhIRg3rx5uHjxIr799ltxfWNtwnbt2oWHHnoI9vb2sLGxgb+/P9544w0At77/fv36AQAmTpwoxl//vQ8dOhS9evVCamoqBg8eDCsrK7HsnW3C6tXW1uKNN96Am5sbrK2t8fjjj+PSpUtq+zTVtvH2Y94rtsbahFVUVODVV18V7yV/f3+89957avc8AMhkMkybNg2bN29Gr169IJfL0bNnT+zYsaPxCicyAnwSRqQH3t7eOHToENLT09GrVy+tjhETEwM3NzcsXLgQhw8fxn/+8x/Y29vj4MGD8PLywrvvvovt27dj5cqV6NWrF1544QWdxL5mzRr07NkTjz/+ODp06IBffvkFr7zyCurq6hAdHQ0AiI+PR0xMDGxsbPDmm28CAFxdXQHceqoyZMgQXL58GVOmTIGXlxcOHjyIuLg45OXlIT4+HgAgCAJGjx6N/fv3Y+rUqejRowc2bdqEyMhInVzHU089hUmTJmHnzp1YsmQJAODo0aM4ePAgxo0bBw8PD1y4cAFr1qzB0KFDcebMGVhZWWHw4MGYPn06PvroI7zxxhvo0aMHAIj/TUhIgI2NDWJjY2FjY4M9e/bg7bffRmlpKVauXNmimCdMmIA33ngDO3fuxOTJkxvd5/Tp03jssccQGBiIRYsWQS6XIzMzEwcOHBDjXLRoEd5++2289NJLGDRoEABgwIAB4jEKCwsRHh6OcePG4fnnnxe/u6YsWbIEMpkMr732GgoKChAfH4/Q0FCkpaWJT+yaozmx3U4QBDz++OPYu3cvJk2ahN69eyMxMRFz5szB5cuX8cEHH6jtv3//fmzcuBGvvPIKbG1t8dFHHyEiIgI5OTlwcnJqdpxEBiMQkc7t3LlTMDU1FUxNTYXg4GBh7ty5QmJiolBdXd1gX29vbyEyMlJcXrdunQBACAsLE+rq6sT1wcHBgkwmE6ZOnSquu3nzpuDh4SEMGTJEXLd3714BgLB3716182RnZwsAhHXr1onr5s+fL9z5a6CysrJBjGFhYUKXLl3U1vXs2VPtvPUWL14sWFtbC3/99Zfa+tdff10wNTUVcnJyBEEQhM2bNwsAhBUrVqhdz6BBgxrE2Zj669ywYUOT+wQFBQkODg53vbZDhw4JAISvv/5aXLdhw4ZG67CpY0yZMkWwsrISqqqq7hpz/Xd79OjRJvdRKBTC/fffLy7f+R198MEHAgDh6tWrTR7j6NGjTdbhkCFDBADC2rVrG93W2L3UqVMnobS0VFz/008/CQCEDz/8UFx3533c1DHvFltkZKTg7e0tLtffI++8847afk899ZQgk8mEzMxMcR0AwdzcXG3diRMnBADCxx9/3OBcRMaAryOJ9ODhhx/GoUOH8Pjjj+PEiRNYsWIFwsLC0KlTJ2zdurVZx5g0aZLaa6j+/ftDEARMmjRJXGdqaoq+ffvif//7n85iv/3JRklJCa5du4YhQ4bgf//7H0pKSu5ZfsOGDRg0aBAcHBxw7do18RMaGora2lrs27cPALB9+3Z06NABL7/8str13N5poaVsbGzUehzefm01NTUoLCyEn58f7O3tcfz48WYd8/ZjlJWV4dq1axg0aBAqKytx7tw5ncd8p/pX3Fu2bNG6EbtcLsfEiRObvf8LL7wAW1tbcfmpp56Cu7s7tm/frtX5m2v79u0wNTXF9OnT1da/+uqrEAQBv/32m9r60NBQ+Pr6isuBgYGws7PT6b8PIl1iEkakJ/369cPGjRtRVFSEP/74A3FxcSgrK8NTTz2FM2fO3LO8l5eX2rJCoQAAeHp6NlhfVFSks7gPHDiA0NBQWFtbw97eHi4uLmKboeYkYefPn8eOHTvg4uKi9gkNDQXwT8eEixcvwt3dHTY2Nmrl/f39dXYt5eXlasnDjRs38Pbbb4vti5ydneHi4oLi4uJmXRtw63XgE088AYVCATs7O7i4uOD5558H0Lz60TTmO40dOxYDBw7Ev//9b7i6umLcuHH46aefNErIOnXqpFEj/K5du6oty2Qy+Pn56XwIjztdvHgRSqWyQX3Uvxq+ePGi2vo7/80AgIODg07/fRDpEtuEEemZubk5+vXrh379+qFbt26YOHEiNmzYgPnz59+1nKmpabPXC7c1Um5qYM/mNBzPysrC8OHD0b17d6xatQqenp4wNzfH9u3b8cEHHzTrD31dXR0efvhhzJ07t9Ht3bp1u+cxdKGmpgZ//fWXWpu8mJgYrFu3DjNnzkRwcDAUCgVkMhnGjRvXrGsrLi7GkCFDYGdnh0WLFsHX1xcWFhY4fvw4XnvttRYPr/D333+jpKQEfn5+Te5jaWmJffv2Ye/evfj111+xY8cO/PjjjwgJCcHOnTubvG/uPIau3e2+a05MutDUeYQ7GvETGQsmYUQGVD8sQV5ent7O4eDgAAANehne+dSgMb/88gtUKhW2bt2q9lRh7969DfZt6o+ur68vysvLxSdfTfH29kZSUhLKy8vVnoZlZGTcM87m+Pnnn3Hjxg2EhYWprYuMjMT7778vrquqqmpQV01dW3JyMgoLC7Fx40YMHjxYXJ+dna2TmL/55hsAUIu5MSYmJhg+fDiGDx+OVatW4d1338Wbb76JvXv3IjQ0VOcj7J8/f15tWRAEZGZmqo1n5uDg0GjP1osXL6JLly7isiaxeXt7Y/fu3SgrK1N7Glb/2tfb27vZxyIyRnwdSaQHe/fubfT/vuvb0OjyldudvL29YWpqKra9qvfpp5/es2z9k4TbYy8pKcG6desa7Gttbd3oH91nnnkGhw4dQmJiYoNtxcXFuHnzJgDgkUcewc2bN9WGv6itrdXJjAInTpzAzJkz4eDgIPboBG5d353fy8cff9zgKaG1tbUY7+0aq5/q6upm1e297NmzB4sXL4aPjw+ee+65Jve7fv16g3X1g56qVCoATcevra+//lqtndrPP/+MvLw8hIeHi+t8fX1x+PBhVFdXi+u2bdvWYCgLTWJ75JFHUFtbqzbMCAB88MEHkMlkaucnao34JIxID2JiYlBZWYknnngC3bt3R3V1NQ4ePIgff/wRnTt31qhRtKYUCgWefvppfPzxx5DJZPD19cW2bdsaDBLbmBEjRsDc3ByjRo3ClClTUF5ejs8++wwdO3Zs8PSuT58+WLNmDd555x34+fmhY8eOCAkJwZw5c7B161Y89thjiIqKQp8+fVBRUYFTp07h559/xoULF+Ds7IxRo0Zh4MCBeP3113HhwgUEBARg48aNGrer+v3331FVVYXa2loUFhbiwIED2Lp1KxQKBTZt2gQ3Nzdx38ceewzffPMNFAoFAgICcOjQIezevbvB8AW9e/eGqakpli9fjpKSEsjlcoSEhGDAgAFwcHBAZGQkpk+fDplMhm+++Ubj112//fYbzp07h5s3b+LKlSvYs2cPdu3aBW9vb2zduhUWFhZNll20aBH27duHRx99FN7e3igoKMCnn34KDw8PPPTQQwBuJUT29vZYu3YtbG1tYW1tjf79+8PHx0ejOOs5OjrioYcewsSJE3HlyhXEx8fDz89PbRiNf//73/j5558xcuRIPPPMM8jKysK3336r1lBe09hGjRqFYcOG4c0338SFCxcQFBSEnTt3YsuWLZg5c2aDYxO1OlJ1yyRqy3777TfhxRdfFLp37y7Y2NgI5ubmgp+fnxATEyNcuXJFbd+mhqi4cxiD+qEK7hyaIDIyUrC2tlZbd/XqVSEiIkKwsrISHBwchClTpgjp6enNGqJi69atQmBgoGBhYSF07txZWL58ufDll18KAITs7Gxxv/z8fOHRRx8VbG1tBQBqwxCUlZUJcXFxgp+fn2Bubi44OzsLAwYMEN577z21YToKCwuFCRMmCHZ2doJCoRAmTJgg/PnnnxoNUVH/MTMzE1xcXITBgwcLS5YsEQoKChqUKSoqEiZOnCg4OzsLNjY2QlhYmHDu3LlGh1f47LPPhC5dugimpqZqw1UcOHBAePDBBwVLS0tBqVSKw4+giSEtblf/3dZ/zM3NBTc3N+Hhhx8WPvzwQ7VhIOrd+R0lJSUJo0ePFpRKpWBubi4olUph/PjxDYYE2bJlixAQECB06NBBrT6HDBki9OzZs9H4mhqi4vvvvxfi4uKEjh07CpaWlsKjjz4qXLx4sUH5999/X+jUqZMgl8uFgQMHCseOHWtwzLvFducQFYJw616aNWuWoFQqBTMzM6Fr167CypUr1YZvEYRbQ1RER0c3iKmpoTOIjIFMENhikYiIiMjQ2CaMiIiISAJMwoiIiIgkwCSMiIiISAJMwoiIiIgkwCSMiIiISAJMwoiIiIgk0OYHa62rq0Nubi5sbW11PpUHERER0e0EQUBZWRmUSiVMTO7+rKvNJ2G5ubnw9PSUOgwiIiJqRy5dugQPD4+77tPmk7D6SV8vXboEOzs7iaMhIiKitqy0tBSenp5qk843pc0nYfWvIO3s7JiEERERkUE0pwkUG+YTERERSYBJGBEREZEEmIQRERERSaDNtwkjIiIyRoIg4ObNm6itrZU6FNKQmZkZTE1NW3wcJmFEREQGVl1djby8PFRWVkodCmlBJpPBw8MDNjY2LToOkzAiIiIDqqurQ3Z2NkxNTaFUKmFubs7BxFsRQRBw9epV/P333+jatWuLnogxCSMiIjKg6upq1NXVwdPTE1ZWVlKHQ1pwcXHBhQsXUFNTwySMqLlyc3NRVFSkcTkHBwcolUo9RERE7dW9prQh46WrJ5dMwqjdyM3NRXd/f5SVl2tc1tbGBucyMpiIERGRzjAJo3ajqKgIZeXl2Lz4efgpnZpdLjO3EGPmfYuioiImYUREpDNMwqjd8VM6oaePq9RhEBFJKioqCsXFxdi8ebPUobRbfCFNRERkRKKioiCTySCTyWBubg4/Pz8sWrQIN2/elDq0u0pISIC9vb3G5TIzM/Hiiy/Cy8sLcrkcnTp1wvDhw/Hdd98Z/TW3FJ+EERERGZmRI0di3bp1UKlU2L59O6Kjo2FmZoa4uDi1/aqrq2Fubi5RlC33xx9/IDQ0FD179sTq1avRvXt3AMCxY8ewevVq9OrVC0FBQY2WrampgZmZmSHD1Tk+CSMiIjIycrkcbm5u8Pb2xssvv4zQ0FBs3boVUVFRGDNmDJYsWQKlUgl/f38AwKlTpxASEgJLS0s4OTnhpZdeQvltnZBqa2sRGxsLe3t7ODk5Ye7cuRAEQe2cnTt3Rnx8vNq63r17Y8GCBeJycXExpkyZAldXV1hYWKBXr17Ytm0bkpOTMXHiRJSUlIhP8W4v1xhBEBAVFYVu3brhwIEDGDVqFLp27YquXbti/Pjx2L9/PwIDAwEAFy5cgEwmw48//oghQ4bAwsIC3333Herq6rBo0SJ4eHhALpejd+/e2LFjh3iO5ORkyGQyFBcXi+vS0tIgk8lw4cIFAP88wdu8eTO6du0KCwsLhIWF4dKlS838trTHJIyIiMjIWVpaorq6GgCQlJSEjIwM7Nq1C9u2bUNFRQXCwsLg4OCAo0ePYsOGDdi9ezemTZsmln///feRkJCAL7/8Evv378f169exadMmjWKoq6tDeHg4Dhw4gG+//RZnzpzBsmXLYGpqigEDBiA+Ph52dnbIy8tDXl4eZs+efdfjpaWl4ezZs5g9e3aTw3XcORTE66+/jhkzZuDs2bMICwvDhx9+iPfffx/vvfceTp48ibCwMDz++OM4f/68RtdWWVmJJUuW4Ouvv8aBAwdQXFyMcePGaXQMbfB1JBERkZESBAFJSUlITExETEwMrl69Cmtra3z++efia8jPPvsMVVVV+Prrr2FtbQ0A+OSTTzBq1CgsX74crq6uiI+PR1xcHJ588kkAwNq1a5GYmKhRLLt378Yff/yBs2fPolu3bgCALl26iNsVCgVkMhnc3Nyadby//voLAMSneQBQUFCgdswVK1bglVdeEZdnzpwpXgMAvPfee3jttdfEhGn58uXYu3cv4uPjsXr16mZfW01NDT755BP0798fAPDVV1+hR48e+OOPP/DAAw80+ziakvRJ2Jo1axAYGAg7OzvY2dkhODgYv/32m7h96NCh4mPN+s/UqVMljJiIiEj/tm3bBhsbG1hYWCA8PBxjx44VX+/dd999au3Azp49i6CgIDEBA4CBAweirq4OGRkZKCkpQV5enphgAECHDh3Qt29fjWJKS0uDh4eHmIDpg5OTE9LS0pCWlgZ7e3vx6V+922MuLS1Fbm4uBg4cqLbPwIEDcfbsWY3O26FDB/Tr109c7t69O+zt7TU+jqYkfRLm4eGBZcuWoWvXrhAEAV999RVGjx6NP//8Ez179gQATJ48GYsWLRLLcIoHIiJq64YNG4Y1a9bA3NwcSqUSHTr88+f69mRLl0xMTBq0E6upqRF/trS01On5unbtCgDIyMjA/fffDwAwNTWFn58fAKhdcz1Nr73+Neft13X7NUlN0idho0aNwiOPPIKuXbuiW7duWLJkCWxsbHD48GFxHysrK7i5uYkfOzs7CSMmIiLSP2tra/j5+cHLy6vRZOR2PXr0wIkTJ1BRUSGuO3DgAExMTODv7w+FQgF3d3ccOXJE3H7z5k2kpqaqHcfFxQV5eXnicmlpKbKzs8XlwMBA/P333+JrxDuZm5ujtra22dd4//33o3v37njvvfdQV1fX7HL17OzsoFQqceDAAbX1Bw4cQEBAAIBb1wRA7brS0tIaHOvmzZs4duyYuJyRkYHi4mL06NFD47g0YTQN82tra/HDDz+goqICwcHB4vrvvvsOzs7O6NWrF+Li4lBZWXnX46hUKpSWlqp9iIiI2qrnnnsOFhYWiIyMRHp6Ovbu3YuYmBhMmDABrq63BqaeMWMGli1bhs2bN+PcuXN45ZVX1HoMAkBISAi++eYb/P777zh16hQiIyPVJqceMmQIBg8ejIiICOzatQvZ2dn47bffxN6InTt3Rnl5OZKSknDt2rV7/r2WyWRYt24dMjIyMHDgQGzduhXnz5/HmTNnsHbtWly9evWek2PPmTMHy5cvx48//oiMjAy8/vrrSEtLw4wZMwAAfn5+8PT0xIIFC3D+/Hn8+uuveP/99xscx8zMDDExMThy5AhSU1MRFRWFBx98UK/twQAjaJh/6tQpBAcHo6qqCjY2Nti0aZOYwT777LPw9vaGUqnEyZMn8dprryEjIwMbN25s8nhLly7FwoULDRU+ERGRpKysrJCYmIgZM2agX79+sLKyQkREBFatWiXu8+qrryIvLw+RkZEwMTHBiy++iCeeeAIlJSXiPnFxccjOzsZjjz0GhUKBxYsXqz0JA4D//ve/mD17NsaPH4+Kigr4+flh2bJlAIABAwZg6tSpGDt2LAoLCzF//vx7DlPx4IMPIjU1Fe+++y6io6ORn58Pa2trBAUF4YMPPsCLL7541/LTp09HSUkJXn31VRQUFCAgIABbt24VX3WamZnh+++/x8svv4zAwED069cP77zzDp5++ukGdfjaa6/h2WefxeXLlzFo0CB88cUX96z7lpIJd74ANrDq6mrk5OSgpKQEP//8Mz7//HOkpKSIidjt9uzZg+HDhyMzMxO+vr6NHk+lUkGlUonLpaWl8PT0RElJCV9ltnOnT59Gr169kP7FDI2mLTqdfQW9Jn2I9PR0sa0iEZG2qqqqkJ2dDR8fH1hYWEgdTruXkJCAmTNnNngyeDd3+w5LS0uhUCialXdI/iSsfkoGAOjTpw+OHj2KDz/8EP/3f//XYN/6nh13S8Lkcjnkcrn+AiYiIiLSAaNpE1avrq5O7UnW7eob07m7uxswIiIiItLU77//DhsbmyY/JPGTsLi4OISHh8PLywtlZWVYv349kpOTkZiYiKysLKxfvx6PPPIInJyccPLkScyaNQuDBw8WpzEgIiIi49S3b99GeyIam6ioKERFRUlybkmTsIKCArzwwgvIy8uDQqFAYGAgEhMT8fDDD+PSpUvYvXs34uPjUVFRAU9PT0REROCtt96SMmQiIiJqBktLS7G5ETVO0iTsbj0PPD09kZKSYsBoiIiIiAzH6NqEEREREbUHTMKIiIiIJMAkjIiIiEgCTMKIiIiIJCD5YK1ERERkeDk5Obh27ZrBzufs7AwvLy+Dna81YBJGRETUzuTk5KB7jx64cY9JtnXJ0soK586e1TgRW716NVauXIn8/HwEBQXh448/1vvE2obCJIyIiKiduXbtGm5UVuK511bC1avxaQB16UpOFr5bPgfXrl3TKAn78ccfERsbi7Vr16J///6Ij49HWFgYMjIy0LFjRz1GbBhMwoiIiNopVy9feHTtKXUYTVq1ahUmT56MiRMnAgDWrl2LX3/9FV9++SVef/11iaNrOTbMJyIiIqNTXV2N1NRUhIaGiutMTEwQGhqKQ4cOSRiZ7jAJIyIiIqNz7do11NbWwtXVVW29q6sr8vPzJYpKt5iEEREREUmASRgREREZHWdnZ5iamuLKlStq669cuQI3NzeJotItJmFERERkdMzNzdGnTx8kJSWJ6+rq6pCUlITg4GAJI9Md9o4kIiJqp67kZBn1eWJjYxEZGYm+ffvigQceQHx8PCoqKsTekq0dkzAiIqJ2xtnZGZZWVvhu+RyDndPSygrOzs4alRk7diyuXr2Kt99+G/n5+ejduzd27NjRoLF+a8UkjIiIqJ3x8vLCubNnW8W0RdOmTcO0adP0EJH0mIQRERG1Q15eXpzLUWJsmE9EREQkAUmTsDVr1iAwMBB2dnaws7NDcHAwfvvtN3F7VVUVoqOj4eTkBBsbG0RERDToqkpERETUGkmahHl4eGDZsmVITU3FsWPHEBISgtGjR+P06dMAgFmzZuGXX37Bhg0bkJKSgtzcXDz55JNShkxERESkE5K2CRs1apTa8pIlS7BmzRocPnwYHh4e+OKLL7B+/XqEhIQAANatW4cePXrg8OHDePDBB6UImYiIiEgnjKZNWG1tLX744QdUVFQgODgYqampqKmpUZu4s3v37vDy8rrrxJ0qlQqlpaVqHyIiIiJjI3kSdurUKdjY2EAul2Pq1KnYtGkTAgICkJ+fD3Nzc9jb26vtf6+JO5cuXQqFQiF+PD099XwFRERERJqTPAnz9/dHWloajhw5gpdffhmRkZE4c+aM1seLi4tDSUmJ+Ll06ZIOoyUiIiLSDcnHCTM3N4efnx8AoE+fPjh69Cg+/PBDjB07FtXV1SguLlZ7GnaviTvlcjnkcrm+wyYiIiJqEcmTsDvV1dVBpVKhT58+MDMzQ1JSEiIiIgAAGRkZyMnJaTMTdxIREUklJyenVYyY35ZJmoTFxcUhPDwcXl5eKCsrw/r165GcnIzExEQoFApMmjQJsbGxcHR0hJ2dHWJiYhAcHMyekURERC2Qk5ODHj26o7LyhsHOaWVlibNnzzU7Edu3bx9WrlyJ1NRU5OXlYdOmTRgzZox+gzQwSZOwgoICvPDCC8jLy4NCoUBgYCASExPx8MMPAwA++OADmJiYICIiAiqVCmFhYfj000+lDJmIiKjVu3btGiorb+DbN55BDy8XvZ/vbM5VPP/uT7h27Vqzk7CKigoEBQXhxRdfbLNjhEqahH3xxRd33W5hYYHVq1dj9erVBoqIiIio/ejh5YJ/deskdRiNCg8PR3h4uNRh6JXkvSOJiIiI2iMmYUREREQSYBJGREREJAEmYUREREQSYBJGREREJAGjG6yViIiIDONszlWjPU95eTkyMzPF5ezsbKSlpcHR0bHNDPrKJIyIiKidcXZ2hpWVJZ5/9yeDndPKyhLOzs7N3v/YsWMYNmyYuBwbGwsAiIyMREJCgq7DkwSTMCIionbGy8sLZ8+eM+ppi4YOHQpBEPQYkfSYhBEREbVDXl5ebea1XmvFhvlEREREEmASRkRERCQBJmFEREREEmASRkREJIG23ui8LdPVd8ckjIiIyIDMzMwAAJWVlRJHQtqqrq4GAJiamrboOOwdSUREZECmpqawt7dHQUEBAMDKygoymUziqKi56urqcPXqVVhZWaFDh5alUUzCiIiIDMzNzQ0AxESMWhcTExN4eXm1OHlmEkatVm5uLoqKipq9/+3TX2hD0/IqlQpyuVyjMg4ODlAqlRqVIaLWRyaTwd3dHR07dkRNTY3U4ZCGzM3NYWLS8hZdTMKoVcrNzUV3f3+UlZdrXLaiskKj/QuKyiEDMGbMGI3KmciAOg3bbtra2OBcRgYTMaJ2wtTUtMXtiqj1YhJGrVJRURHKysuxefHz8FM6NavMnj+zMP2TbahSqTQ6V2mlCgKAhNmj0LdHF43OpUmZzNxCjJn3LYqKipiEERG1A5ImYUuXLsXGjRtx7tw5WFpaYsCAAVi+fDn8/f3FfYYOHYqUlBS1clOmTMHatWsNHS4ZIT+lE3r6uDZr38zLhS06l4+bvcbn0qQMERG1L5IOUZGSkoLo6GgcPnwYu3btQk1NDUaMGIGKCvXXRZMnT0ZeXp74WbFihUQRExEREemGpE/CduzYobackJCAjh07IjU1FYMHDxbXW1lZiT1JiIiIiNoCo2oTVlJSAgBwdHRUW//dd9/h22+/hZubG0aNGoV58+bBysqq0WOoVCqobmvzU1paqr+AqQFNeyzWY69AIiJqb4wmCaurq8PMmTMxcOBA9OrVS1z/7LPPwtvbG0qlEidPnsRrr72GjIwMbNy4sdHjLF26FAsXLjRU2HSblvRYZK9AIiJqb4wmCYuOjkZ6ejr279+vtv6ll14Sf77vvvvg7u6O4cOHIysrC76+vg2OExcXh9jYWHG5tLQUnp6e+gucRNr0WATYK5CIiNono0jCpk2bhm3btmHfvn3w8PC46779+/cHcGvgzMaSMLlcrvEAmaRbmvRYJCIiaq8kTcIEQUBMTAw2bdqE5ORk+Pj43LNMWloaAMDd3V3P0RERERHpj6RJWHR0NNavX48tW7bA1tYW+fn5AACFQgFLS0tkZWVh/fr1eOSRR+Dk5ISTJ09i1qxZGDx4MAIDA6UMnYiIiKhFJE3C1qxZA+DWgKy3W7duHaKiomBubo7du3cjPj4eFRUV8PT0REREBN566y0JoiUiIiLSHclfR96Np6dng9HyiYiIiNoCSUfMJyIiImqvmIQRERERSYBJGBEREZEEmIQRERERSYBJGBEREZEEmIQRERERSYBJGBEREZEEmIQRERERSUCrJKxLly4oLCxssL64uBhdunRpcVBEREREbZ1WSdiFCxdQW1vbYL1KpcLly5dbHBQRERFRW6fRtEVbt24Vf05MTIRCoRCXa2trkZSUhM6dO+ssOCIiIqK2SqMkbMyYMQAAmUyGyMhItW1mZmbo3Lkz3n//fZ0FR0RERNRWaZSE1dXVAQB8fHxw9OhRODs76yUoIiIiorZOoySsXnZ2tq7jICIiImpXtErCACApKQlJSUkoKCgQn5DV+/LLL1scGBEREVFbplUStnDhQixatAh9+/aFu7s7ZDKZruMiIiIiatO0SsLWrl2LhIQETJgwQdfxEBEREbULWo0TVl1djQEDBug6FiIiIqJ2Q6sk7N///jfWr1+v61iIiIiI2g2tXkdWVVXhP//5D3bv3o3AwECYmZmpbV+1alWzjrN06VJs3LgR586dg6WlJQYMGIDly5fD399f7VyvvvoqfvjhB6hUKoSFheHTTz+Fq6urNqETERERGQWtkrCTJ0+id+/eAID09HS1bZo00k9JSUF0dDT69euHmzdv4o033sCIESNw5swZWFtbAwBmzZqFX3/9FRs2bIBCocC0adPw5JNP4sCBA9qETkRERGQUtErC9u7dq5OT79ixQ205ISEBHTt2RGpqKgYPHoySkhJ88cUXWL9+PUJCQgAA69atQ48ePXD48GE8+OCDOomDiIiIyNC0HidMH0pKSgAAjo6OAIDU1FTU1NQgNDRU3Kd79+7w8vLCoUOHGk3CVCoVVCqVuFxaWqrnqElXMjMz9bIvERGRMdIqCRs2bNhdXzvu2bNH42PW1dVh5syZGDhwIHr16gUAyM/Ph7m5Oezt7dX2dXV1RX5+fqPHWbp0KRYuXKjx+Uk6BUXlkOGfuUk1UVFZofN4iIiIDEGrJKy+PVi9mpoapKWlIT09vcHE3s0VHR2N9PR07N+/X6vy9eLi4hAbGysul5aWwtPTs0XHJP0qrVRBAJAwexT69ujSrDJ7/szC9E+2oeq2p55EREStiVZJ2AcffNDo+gULFqC8vFzj402bNg3btm3Dvn374OHhIa53c3NDdXU1iouL1Z6GXblyBW5ubo0eSy6XQy6XaxwDSc/HzR49fZrX6zXzcqGeoyEiItIvrcYJa8rzzz+v0byRgiBg2rRp2LRpE/bs2QMfHx+17X369IGZmRmSkpLEdRkZGcjJyUFwcLDO4iYiIiIyNJ02zD906BAsLCyavX90dDTWr1+PLVu2wNbWVmznpVAoYGlpCYVCgUmTJiE2NhaOjo6ws7NDTEwMgoOD2TOSiIiIWjWtkrAnn3xSbVkQBOTl5eHYsWOYN29es4+zZs0aAMDQoUPV1q9btw5RUVEAbr36NDExQUREhNpgrUREREStmVZJmEKhUFs2MTGBv78/Fi1ahBEjRjT7OIIg3HMfCwsLrF69GqtXr9Y4TiIiIiJjpVUStm7dOl3HQURERNSutKhNWGpqKs6ePQsA6NmzJ+6//36dBEVERETU1mmVhBUUFGDcuHFITk4Wh44oLi7GsGHD8MMPP8DFxUWXMRIRERG1OVoNURETE4OysjKcPn0a169fx/Xr15Geno7S0lJMnz5d1zESERERtTlaPQnbsWMHdu/ejR49eojrAgICsHr1ao0a5hMRERG1V1o9Caurq4OZmVmD9WZmZqirq2txUERERERtnVZJWEhICGbMmIHc3Fxx3eXLlzFr1iwMHz5cZ8ERERERtVVaJWGffPIJSktL0blzZ/j6+sLX1xc+Pj4oLS3Fxx9/rOsYiYiIiNocrdqEeXp64vjx49i9ezfOnTsHAOjRowdCQ0N1GhwRERFRW6XRk7A9e/YgICAApaWlkMlkePjhhxETE4OYmBj069cPPXv2xO+//66vWImIiIjaDI2SsPj4eEyePBl2dnYNtikUCkyZMgWrVq3SWXBEREREbZVGSdiJEycwcuTIJrePGDECqampLQ6KiIiIqK3TKAm7cuVKo0NT1OvQoQOuXr3a4qCIiIiI2jqNkrBOnTohPT29ye0nT56Eu7t7i4MiIiIiaus0SsIeeeQRzJs3D1VVVQ223bhxA/Pnz8djjz2ms+CIiIiI2iqNhqh46623sHHjRnTr1g3Tpk2Dv78/AODcuXNYvXo1amtr8eabb+olUCIiIqK2RKMkzNXVFQcPHsTLL7+MuLg4CIIAAJDJZAgLC8Pq1avh6uqql0CJ2ovMzEyNyzg4OECpVOohGiIi0heNB2v19vbG9u3bUVRUhMzMTAiCgK5du8LBwUEf8RG1GwVF5ZABGDNmjMZlbW1scC4jg4kYEVErotWI+cCt//Pu16+fLmMhatdKK1UQACTMHoW+Pbo0u1xmbiHGzPsWRUVFTMKIiFoRrZMwXdi3bx9WrlyJ1NRU5OXlYdOmTWpPAaKiovDVV1+plQkLC8OOHTsMHCmR4fi42aOnD1/rExG1dVpN4K0rFRUVCAoKwurVq5vcZ+TIkcjLyxM/33//vQEjJCIiItIPSZ+EhYeHIzw8/K77yOVyuLm5GSgiIiIiIsOQNAlrjuTkZHTs2BEODg4ICQnBO++8Aycnpyb3V6lUUKlU4nJpaakhwiQiI5Cbm4uioiKNymjbs1Sbc6lUKsjlco3KsOcrUdtl1EnYyJEj8eSTT8LHxwdZWVl44403EB4ejkOHDsHU1LTRMkuXLsXChQsNHCkRSS03Nxf+/t1RXl6mUTkbG1tkZJzTKNHR9lyQyYD/P7SPPuMjotbBqJOwcePGiT/fd999CAwMhK+vL5KTkzF8+PBGy8TFxSE2NlZcLi0thaenp95jJSJpFRUVoby8DC8uWA1npXezylzLvYgvF0Rr3LNUm3OdTzuETZ8uQcSsd+HbI1Cv8RFR62DUSdidunTpAmdnZ2RmZjaZhMnlco0f9xNR2+Gs9IZb565Gd65ruRcBAA5uHgaLj4iMm6S9IzX1999/o7CwkJOEExERUasn6ZOw8vJytSlasrOzkZaWBkdHRzg6OmLhwoWIiIiAm5sbsrKyMHfuXPj5+SEsLEzCqFsnQzQi1ma6HWofDNlgnoiotZA0CTt27BiGDRsmLte35YqMjMSaNWtw8uRJfPXVVyguLoZSqcSIESOwePFivm7UUG5uLrr7+6OsvFyjciYyoE6zNsQAgIrKCs0LUZtlyAbzREStiaRJ2NChQ8VJwBuTmJhowGjarqKiIpSVl2Pz4ufhp2x6eI/b7fkzC9M/2abRFDr1ZapuGyKEyJAN5omIWpNW1TCfWsZP6dTs6XAyLxcC0GwKnfoyRI0xZIN5IqLWoFU1zCciIiJqK5iEEREREUmAryOJyGhp0uO2LffO1eba2LuUyPgxCSMio1NWXAhAhjFjxmhctqKi7fTObUk9sHcpkfFjEkZERkdVWQ5A0GiKn/ppgVTVbad3rjb1ALB3KVFrwSSMiIyWJlP81E8L1BZxqiOitokN84mIiIgkwCSMiIiISAJ8HdnKaDMHX1vuNUZETdP03z57VBqeNr/TAX5XbQWTsFZE2zkg63FOR6L2QdtelexRaVjazqsK8LtqK5iEtSLazAEJcE5HovZGm16V7FFpeNrMqwrwu2pLmIS1QprMAQlwTkei9oq9KlsHzqvafrFhPhEREZEEmIQRERERSYCvI4naKW16ZbFHVtvHHpVEhsMkjKgd0ranra2NDc5lZPCPbhvEHpVEhsckjKgd0qanbWZuIcbM+5Y9stoo9qgkMjxJk7B9+/Zh5cqVSE1NRV5eHjZt2qT2f2GCIGD+/Pn47LPPUFxcjIEDB2LNmjXo2pW9SIh0QdOettT2sUclkeFI2jC/oqICQUFBWL16daPbV6xYgY8++ghr167FkSNHYG1tjbCwMFRVVRk4UiIiIiLdkvRJWHh4OMLDwxvdJggC4uPj8dZbb2H06NEAgK+//hqurq7YvHkzxo0bZ8hQiYiIiHTKaNuEZWdnIz8/H6GhoeI6hUKB/v3749ChQ00mYSqVCqrbRoYvLS3Ve6xExmDPnj3N7tmWk5Oj52iIiOhejDYJy8/PBwC4uqq3V3F1dRW3NWbp0qVYuHChXmMjMiZnLlwBIMP06dM1LCnDhbwCtgkjIpKI0SZh2oqLi0NsbKy4XFpaCk9PTwkjItKv/KJbvdqGjRkHX5/OzSpzLusC9m/9AdeKtZsMnoiIWs5okzA3NzcAwJUrV+Du7i6uv3LlCnr37t1kOblcDrlcru/wiIxOp07uCOjRvF5tJZWczJ2ISGpGO22Rj48P3NzckJSUJK4rLS3FkSNHEBwcLGFkRERERC0n6ZOw8vJytYbE2dnZSEtLg6OjI7y8vDBz5ky888476Nq1K3x8fDBv3jwolUqNR3QmIiIiMjaSJmHHjh3DsGHDxOX6tlyRkZFISEjA3LlzUVFRgZdeegnFxcV46KGHsGPHDlhYWEgVMhER3UHT+SaBtjnnpKbzsWpTb9S2SJqEDR06FIIgNLldJpNh0aJFWLRokQGjIiKi5tB2vkmg7c05mZubC3//7igvL9O4bEVFhR4iotbAaBvmExGRcdNmvkmgbc45WVRUhPLyMry4YDWcld7NKnM+7RA2fboEqmp2lGmvmIQREVGLcL7JfzgrvZtdF9dyL+o5GjJ2Rts7koiIiKgtYxJGREREJAG+jqRWLTu/GE7ZV5q1b05BscHPZcj4SHua9lJjrzZpaNr7EGibvTCp7WASRq3S9dJKADJEvfeLhiVluFZSabBzGSI+0l5LevcB7NVmSNr2PmxrvTCpbWESRq1SeVU1AAGPjZ+AAP/mNYI9cfo8Ejd8g7Ib1Xo/14GDR3Fg5y8YNmYc+vUO0Gt8pD1te/exV5vhadP7sC32wqS2hUkYtWr2js5wU3Zq1r7ZeZq9xmjJuawV5wEAtg6OBouPtKdp7z72apOOJr0PiYwdG+YTERERSYBJGBEREZEE+DpSQpxnTBqXr5XhdDN7LAJAXqHm05C0hLHHR6QrmvxOM/TvP/5+JkNgEiaR3NxcdPf3R1l5ucZlKyrZI0sbNyrKAcjwZkIy3kxI1ry8qkbnMakd38jjI9KVlvRKNUSPVM4DSYbCJEwiRUVFKCsvx+bFz8NP6dSsMnv+zML0T7ahSsUeWdqoVlUBEDTqsQj809OxprZOf8HB+OMj0hVteqUaskcq54EkQ2ESJjE/pRN6+rg2a9/My4V6jqZ90KTHIvBPT0dDMfb4iHRFk16pUvRI5TyQpG9smE9EREQkASZhRERERBLg60gi0gjnWSQi0g0mYUTULAVFt3pvcp5FIiLdYBJGRM1SWqkC51kkItIdo07CFixYgIULF6qt8/f3x7lz5ySKiIg4zyIRkW4YdRIGAD179sTu3bvF5Q4djD5kIiIionsy+oymQ4cOcHNzkzoMIiIiIp0y+iTs/PnzUCqVsLCwQHBwMJYuXQovL68m91epVFDdNqJ8aWmpIcIkapU0macyp6BYv8EQ6Ql79JKxMuokrH///khISIC/vz/y8vKwcOFCDBo0COnp6bC1tW20zNKlSxu0IyMidS2Zp/LGjUq9xESkay2ZoxJgj17SP6NOwsLDw8WfAwMD0b9/f3h7e+Onn37CpEmTGi0TFxeH2NhYcbm0tBSenp56j5WoNdFmnso/j5/Arq3/RXX1Tf0GR6Qj2sxRCbBHLxmOUSdhd7K3t0e3bt3u+qhYLpdDLpcbMCqi1kuTeSrtstjLkVon9uglY9Wqpi0qLy9HVlYW3N3dpQ6FiIiIqEWMOgmbPXs2UlJScOHCBRw8eBBPPPEETE1NMX78eKlDIyIiImoRo34d+ffff2P8+PEoLCyEi4sLHnroIRw+fBguLi5Sh0ZERNSq5ObmoqioSKMyDg4OUCqVeoqIjDoJ++GHH6QOgYiIqNXLzc2Fv393lJeXaVTOxsYWGRnnmIjpiVEnYURERNRyRUVFKC8vw4sLVsNZ6d2sMtdyL+LLBdEoKipiEqYnTMKIiIjaCWelt0Y9RUm/jLphPhEREVFbxSSMiIiISAJ8HUlEGikvK8PVqwXN3r+kpESP0RCRMWJPzOZhEkZEzVJVc2u6ot/3/46Dx081u1x14SUAnHOSqL1gT8zmYxJGRM1SU1sHAOjv3wkP3N+r2eVSj5li57nfOeckUTvBnpjNxySMiDRibWkGFwebZu9va2Gux2iIyFixJ+a9sWE+ERERkQSYhBERERFJgK8jdUTTniCZmZl6jKZx2fnFcMq+0qx9cwqK9RtMIzSJL69QswafRESkHU3/Xknx9621YhKmA7m5ueju74+y8nKNy1ZUVughInXXSysByBD13i8alpThWon+e7RpHx9wQ1Wj+4CIiAhlxYUAZBgzZoxW5Ssq9P/3rbVjEqYDRUVFKCsvx+bFz8NP6dSsMnv+zML0T7ahSqXSc3RAeVU1AAGPjZ+AAP/mNZI8cfo8Ejd8g7Ib1foNDtrFd+DgURzY+YvYY4+IiHRLVVkOQEDErHfh2yOw2eXOpx3Cpk+XQFWt/79vrR2TMB3yUzqhp49rs/bNvFyo52gasnd0hpuyU7P2zc7TbJA9XdAkPmvFeT1HQ0REAODg5qFRL8druRf1GE3bwob5RERERBJgEkZEREQkAb6ObIU06UUItKwn4eVrZThtxD0qqfXQZM5JzjdJdG+a9EJsLT0WNY1TpVJBLpdrfB5jmaeSSVgr0pJehIBmPQlvVJQDkOHNhGS8mZCswVkM06OSWg9t5pzkfJNETWtJr0Vj7bGo9TXJZIAgaHw+Y5mnkklYK6JNL0JAu56E1aoqAAKGjRmHfr0DmlXGkD0qqfXQZs5JzjdJ1DRtei0ae4/FllyTpr03jWmeylaRhK1evRorV65Efn4+goKC8PHHH+OBBx6QOizJaNKLEGhZT0JbB0ej7lFJrYcmc05yvkmie9Ok12Jr6bGozTVp2nvTmBh9w/wff/wRsbGxmD9/Po4fP46goCCEhYWhoKB5bUuIiIiIjJHRJ2GrVq3C5MmTMXHiRAQEBGDt2rWwsrLCl19+KXVoRERERFoz6teR1dXVSE1NRVxcnLjOxMQEoaGhOHToUKNlVCoVVLeNQl/fy6q0tFRvcZb//+mKfj1yDiez8ppV5thflwEAu49n4/L15r2jP5Zxq8zFCzkw0aAdYn7erZguX7qMNMvm9SLRpsylnBwAwKEzl2Funtbs+LS5LkNdkyHPxfj+8fflW/fEpczTsJSbNavMpbO3Gv1nnzmB6qobeitjyHMxPsbH+HQfX/GVW79fysvL9ZIb1B9TaE6HAcGIXb58WQAgHDx4UG39nDlzhAceeKDRMvPnzxcA8MMPP/zwww8//Ej2uXTp0j3zHKN+EqaNuLg4xMbGist1dXW4fv06nJycIJPJJIzMcEpLS+Hp6YlLly7Bzs5O6nAkxbr4B+viH6yLf7Au/sG6uIX18A9t6kIQBJSVlTWr56VRJ2HOzs4wNTXFlSvqg4VeuXIFbm5ujZaRy+UNBm6zt7fXV4hGzc7Ort3/A6rHuvgH6+IfrIt/sC7+wbq4hfXwD03rQqFQNGs/o26Yb25ujj59+iApKUlcV1dXh6SkJAQHB0sYGREREVHLGPWTMACIjY1FZGQk+vbtiwceeADx8fGoqKjAxIkTpQ6NiIiISGtGn4SNHTsWV69exdtvv438/Hz07t0bO3bsgKurq9ShGS25XI758+drNZ9WW8O6+Afr4h+si3+wLv7BuriF9fAPfdeFTBC0mHSJiIiIiFrEqNuEEREREbVVTMKIiIiIJMAkjIiIiEgCTMKIiIiIJMAkrA1ZsGABZDKZ2qd79+5Sh2UQ+/btw6hRo6BUKiGTybB582a17YIg4O2334a7uzssLS0RGhqK8+fPSxOsnt2rLqKiohrcJyNHjpQmWD1aunQp+vXrB1tbW3Ts2BFjxoxBRkaG2j5VVVWIjo6Gk5MTbGxsEBER0WBw6LagOXUxdOjQBvfF1KlTJYpYf9asWYPAwEBx8M3g4GD89ttv4vb2ck8A966L9nJP3GnZsmWQyWSYOXOmuE5f9wWTsDamZ8+eyMvLEz/79++XOiSDqKioQFBQEFavXt3o9hUrVuCjjz7C2rVrceTIEVhbWyMsLAxVVVUGjlT/7lUXADBy5Ei1++T77783YISGkZKSgujoaBw+fBi7du1CTU0NRowYgYqKCnGfWbNm4ZdffsGGDRuQkpKC3NxcPPnkkxJGrR/NqQsAmDx5stp9sWLFCoki1h8PDw8sW7YMqampOHbsGEJCQjB69GicPn0aQPu5J4B71wXQPu6J2x09ehT/93//h8DAQLX1ersvWjzLNhmN+fPnC0FBQVKHITkAwqZNm8Tluro6wc3NTVi5cqW4rri4WJDL5cL3338vQYSGc2ddCIIgREZGCqNHj5YkHikVFBQIAISUlBRBEG7dA2ZmZsKGDRvEfc6ePSsAEA4dOiRVmAZxZ10IgiAMGTJEmDFjhnRBScjBwUH4/PPP2/U9Ua++LgSh/d0TZWVlQteuXYVdu3apXbs+7ws+CWtjzp8/D6VSiS5duuC5555DTk6O1CFJLjs7G/n5+QgNDRXXKRQK9O/fH4cOHZIwMukkJyejY8eO8Pf3x8svv4zCwkKpQ9K7kpISAICjoyMAIDU1FTU1NWr3Rffu3eHl5dXm74s766Led999B2dnZ/Tq1QtxcXGorKyUIjyDqa2txQ8//ICKigoEBwe363vizrqo157uiejoaDz66KNq3z+g398VRj9iPjVf//79kZCQAH9/f+Tl5WHhwoUYNGgQ0tPTYWtrK3V4ksnPzweABrMsuLq6itvak5EjR+LJJ5+Ej48PsrKy8MYbbyA8PByHDh2Cqamp1OHpRV1dHWbOnImBAweiV69eAG7dF+bm5rC3t1fbt63fF43VBQA8++yz8Pb2hlKpxMmTJ/Haa68hIyMDGzdulDBa/Th16hSCg4NRVVUFGxsbbNq0CQEBAUhLS2t390RTdQG0r3vihx9+wPHjx3H06NEG2/T5u4JJWBsSHh4u/hwYGIj+/fvD29sbP/30EyZNmiRhZGRMxo0bJ/583333ITAwEL6+vkhOTsbw4cMljEx/oqOjkZ6e3m7aSN5NU3Xx0ksviT/fd999cHd3x/Dhw5GVlQVfX19Dh6lX/v7+SEtLQ0lJCX7++WdERkYiJSVF6rAk0VRdBAQEtJt74tKlS5gxYwZ27doFCwsLg56bryPbMHt7e3Tr1g2ZmZlShyIpNzc3AGjQk+XKlSvitvasS5cucHZ2brP3ybRp07Bt2zbs3bsXHh4e4no3NzdUV1ejuLhYbf+2fF80VReN6d+/PwC0yfvC3Nwcfn5+6NOnD5YuXYqgoCB8+OGH7fKeaKouGtNW74nU1FQUFBTgX//6Fzp06IAOHTogJSUFH330ETp06ABXV1e93RdMwtqw8vJyZGVlwd3dXepQJOXj4wM3NzckJSWJ60pLS3HkyBG1tg/t1d9//43CwsI2d58IgoBp06Zh06ZN2LNnD3x8fNS29+nTB2ZmZmr3RUZGBnJyctrcfXGvumhMWloaALS5+6IxdXV1UKlU7eqeaEp9XTSmrd4Tw4cPx6lTp5CWliZ++vbti+eee078WV/3BV9HtiGzZ8/GqFGj4O3tjdzcXMyfPx+mpqYYP3681KHpXXl5udr/nWVnZyMtLQ2Ojo7w8vLCzJkz8c4776Br167w8fHBvHnzoFQqMWbMGOmC1pO71YWjoyMWLlyIiIgIuLm5ISsrC3PnzoWfnx/CwsIkjFr3oqOjsX79emzZsgW2trZi2w2FQgFLS0soFApMmjQJsbGxcHR0hJ2dHWJiYhAcHIwHH3xQ4uh16151kZWVhfXr1+ORRx6Bk5MTTp48iVmzZmHw4MENuuq3dnFxcQgPD4eXlxfKysqwfv16JCcnIzExsV3dE8Dd66I93RO2trZq7SMBwNraGk5OTuJ6vd0XLepbSUZl7Nixgru7u2Bubi506tRJGDt2rJCZmSl1WAaxd+9eAUCDT2RkpCAIt4apmDdvnuDq6irI5XJh+PDhQkZGhrRB68nd6qKyslIYMWKE4OLiIpiZmQne3t7C5MmThfz8fKnD1rnG6gCAsG7dOnGfGzduCK+88org4OAgWFlZCU888YSQl5cnXdB6cq+6yMnJEQYPHiw4OjoKcrlc8PPzE+bMmSOUlJRIG7gevPjii4K3t7dgbm4uuLi4CMOHDxd27twpbm8v94Qg3L0u2tM90Zg7h+fQ130hEwRBaFkaR0RERESaYpswIiIiIgkwCSMiIiKSAJMwIiIiIgkwCSMiIiKSAJMwIiIiIgkwCSMiIiKSAJMwIiIiIgkwCSMiIiKSAJMwIjJKUVFRBplWqnPnzoiPj9f7eYzZ0KFDMXPmTKnDIGp3mIQRUbNFRUVBJpNBJpPB3Nwcfn5+WLRoEW7evCl1aPeUkJAAe3v7BuuPHj2Kl156Se/n/+yzzxAUFAQbGxvY29vj/vvvx9KlS/V+XiIyXpzAm4g0MnLkSKxbtw4qlQrbt29HdHQ0zMzMEBcX12Df6upqmJubSxBl87m4uOj9HF9++SVmzpyJjz76CEOGDIFKpcLJkyeRnp6u93MTkfHikzAi0ohcLoebmxu8vb3x8ssvIzQ0FFu3bgXwzyvEJUuWQKlUwt/fHwBw6tQphISEwNLSEk5OTnjppZdQXl4uHrO2thaxsbGwt7eHk5MT5s6dizuntW3stWHv3r2xYMECcbm4uBhTpkyBq6srLCws0KtXL2zbtg3JycmYOHEiSkpKxCd59eXuPG5OTg5Gjx4NGxsb2NnZ4ZlnnsGVK1fE7QsWLEDv3r3xzTffoHPnzlAoFBg3bhzKysqarLOtW7fimWeewaRJk+Dn54eePXti/PjxWLJkibhPfd0tXLgQLi4usLOzw9SpU1FdXS3uU1dXh6VLl8LHxweWlpYICgrCzz//rHau9PR0hIeHw8bGBq6urpgwYQKuXbsmbq+oqMALL7wAGxsbuLu74/33328ybiLSLyZhRNQilpaWaolCUlISMjIysGvXLmzbtg0VFRUICwuDg4MDjh49ig0bNmD37t2YNm2aWOb9999HQkICvvzyS+zfvx/Xr1/Hpk2bNIqjrq4O4eHhOHDgAL799lucOXMGy5Ytg6mpKQYMGID4+HjY2dkhLy8PeXl5mD17dqPHGD16NK5fv46UlBTs2rUL//vf/zB27Fi1/bKysrB582Zs27YN27ZtQ0pKCpYtW9ZkbG5ubjh8+DAuXrx412tISkrC2bNnkZycjO+//x4bN27EwoULxe1Lly7F119/jbVr1+L06dOYNWsWnn/+eaSkpAC4lYSGhITg/vvvx7Fjx7Bjxw5cuXIFzzzzjHiMOXPmICUlBVu2bMHOnTuRnJyM48ePN6uOiUjHBCKiZoqMjBRGjx4tCIIg1NXVCbt27RLkcrkwe/Zscburq6ugUqnEMv/5z38EBwcHoby8XFz366+/CiYmJkJ+fr4gCILg7u4urFixQtxeU1MjeHh4iOcSBEHw9vYWPvjgA7V4goKChPnz5wuCIAiJiYmCiYmJkJGR0Wjs69atExQKRYP1tx93586dgqmpqZCTkyNuP336tABA+OOPPwRBEIT58+cLVlZWQmlpqbjPnDlzhP79+zd6XkEQhNzcXOHBBx8UAAjdunUTIiMjhR9//FGora0V94mMjBQcHR2FiooKcd2aNWsEGxsboba2VqiqqhKsrKyEgwcPqh170qRJwvjx4wVBEITFixcLI0aMUNt+6dIlAYCQkZEhlJWVCebm5sJPP/0kbi8sLBQsLS2FGTNmNBk/EekH24QRkUa2bdsGGxsb1NTUoK6uDs8++6zaK8H77rtPrR3Y2bNnERQUBGtra3HdwIEDUVdXh4yMDFhYWCAvLw/9+/cXt3fo0AF9+/Zt8ErybtLS0uDh4YFu3bppfW1nz56Fp6cnPD09xXUBAQGwt7fH2bNn0a9fPwC3XmHa2tqK+7i7u6OgoKDJ47q7u+PQoUNIT0/Hvn37cPDgQURGRuLzzz/Hjh07YGJy66VEUFAQrKysxHLBwcEoLy/HpUuXUF5ejsrKSjz88MNqx66ursb9998PADhx4gT27t0LGxubBjFkZWXhxo0bqK6uVqtrR0dH8bUxERkWkzAi0siwYcOwZs0amJubQ6lUokMH9V8jtydbumRiYtIgKaupqRF/trS01Mt5G2NmZqa2LJPJUFdXd89yvXr1Qq9evfDKK69g6tSpGDRoEFJSUjBs2LB7lq1vQ/frr7+iU6dOatvkcrm4z6hRo7B8+fIG5d3d3ZGZmXnP8xCR4bBNGBFpxNraGn5+fvDy8mqQgDWmR48eOHHiBCoqKsR1Bw4cgImJCfz9/aFQKODu7o4jR46I22/evInU1FS147i4uCAvL09cLi0tRXZ2trgcGBiIv//+G3/99VejcZibm6O2tvaesV66dAmXLl0S1505cwbFxcUICAi457Vqov54t9fLiRMncOPGDXH58OHDsLGxgaenJwICAiCXy5GTkwM/Pz+1T/2Tu3/96184ffo0Onfu3GAfa2tr+Pr6wszMTK2ui4qKmqwzItIvJmFEpFfPPfccLCwsEBkZifT0dOzduxcxMTGYMGECXF1dAQAzZszAsmXLsHnzZpw7dw6vvPIKiouL1Y4TEhKCb775Br///jtOnTqFyMhImJqaituHDBmCwYMHIyIiArt27UJ2djZ+++037NixA8CtV4jl5eVISkrCtWvXUFlZ2SDW0NBQ3HfffXjuuedw/Phx/PHHH3jhhRcwZMgQ9O3bV+s6ePnll7F48WIcOHAAFy9exOHDh/HCCy/AxcUFwcHB4n7V1dWYNGkSzpw5g+3bt2P+/PmYNm0aTExMYGtri9mzZ2PWrFn46quvkJWVhePHj+Pjjz/GV199BQCIjo7G9evXMX78eBw9ehRZWVlITEzExIkTUVtbCxsbG0yaNAlz5szBnj17kJ6ejqioKPF1KBEZFv/lEZFeWVlZITExEdevX0e/fv3w1FNPYfjw4fjkk0/EfV599VVMmDABkZGRCA4Ohq2tLZ544gm148TFxWHIkCF47LHH8Oijj2LMmDHw9fVV2+e///0v+vXrh/HjxyMgIABz584Vn34NGDAAU6dOxdixY+Hi4oIVK1Y0iFUmk2HLli1wcHDA4MGDERoaii5duuDHH39sUR2Ehobi8OHDePrpp9GtWzdERETAwsICSUlJcHJyEvcbPnw4unbtisGDB2Ps2LF4/PHH1drbLV68GPPmzcPSpUvRo0cPjBw5Er/++it8fHwAAEqlEgcOHEBtbS1GjBiB++67DzNnzoS9vb2YaK1cuRKDBg3CqFGjEBoaioceegh9+vRp0fURkXZkgiYtX4mISC+ioqJQXFyMzZs3Sx0KERkIn4QRERERSYBJGBEREZEE+DqSiIiISAJ8EkZEREQkASZhRERERBJgEkZEREQkASZhRERERBJgEkZEREQkASZhRERERBJgEkZEREQkASZhRERERBL4f6vx0sy7x4SWAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 700x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7, 3))\n",
    "sns.histplot(\n",
    "    data=simulated_dataset, \n",
    "    x=\"Production_Speed\", \n",
    "    hue=\"Product_Group\",\n",
    "    binwidth=1\n",
    ")\n",
    "plt.xlabel(\"Production Speed\")\n",
    "plt.title(\"Simulated Data Distribution\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.pairplot(simulated_dataset.iloc[:, 0:], hue=\"Product_Group\", corner=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_dataset[\"Product_Group\"] = (simulated_dataset[\"Product_Group\"]\n",
    "                                      .map({\n",
    "                                          \"Product_Group_1\": 0, \n",
    "                                          \"Product_Group_2\": 1\n",
    "                                          })\n",
    "                                    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameter Masking for Missing Data\n",
    "\n",
    "To handle the missing data in our hierarchical model, we employ a technique called parameter masking. This approach allows us to effectively \"turn off\" certain parameters for specific groups where they are not applicable. Here's how it works in our model:\n",
    "\n",
    "Key aspects of this approach:\n",
    "\n",
    "1. We use an `parameter_mask` array to specify which parameters are relevant for each product group.\n",
    "2. The input data `X` is masked to replace NaNs with zeros, and then multiplied by the indicators.\n",
    "3. The weights `w` are also masked using the indicators, ensuring that irrelevant parameters don't contribute to the predictions.\n",
    "\n",
    "This approach allows the model to learn parameters even when data is missing for some groups, by effectively ignoring those parameters in the relevant calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameter_mask = jnp.array([\n",
    "    [1, 1, 0],\n",
    "    [1, 1, 1]\n",
    "])\n",
    "config_idx = jnp.array(simulated_dataset[\"Product_Group\"].values)\n",
    "X = jnp.array(simulated_dataset.iloc[:, 1:4].values)\n",
    "y = jnp.array(simulated_dataset[\"Production_Speed\"].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(indicators, config_idx, X, y=None):\n",
    "    n_configs, n_params = indicators.shape\n",
    "\n",
    "    # Create a masked version of X where NaNs are treated as zeros based on the indicators\n",
    "    X_masked = jnp.where(jnp.isnan(X), 0.0, X) * indicators[config_idx]\n",
    "\n",
    "    # Group-specific effects\n",
    "    with numpyro.plate(\"config_i\", n_configs):\n",
    "        alpha = numpyro.sample(\"alpha\", dist.Normal(0., 5.))\n",
    "        with numpyro.plate(\"param_i\", n_params):\n",
    "            w = numpyro.sample(\"w\", dist.Normal(0., 5.))\n",
    "\n",
    "    # Compute the weighted sum of features using indicators to zero-out unused parameters\n",
    "    w_masked = jnp.multiply(w.T, indicators)\n",
    "    \n",
    "    eq = alpha[config_idx] + jnp.sum(jnp.multiply(X_masked, w_masked[config_idx]), axis=-1)\n",
    "    mu = numpyro.deterministic(\"mu\", eq)\n",
    "    scale = numpyro.sample(\"scale\", dist.HalfNormal(5.))\n",
    "\n",
    "    with numpyro.plate(\"data\", X.shape[0]):\n",
    "        numpyro.sample(\"obs\", dist.Normal(mu, scale), obs=y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 12.1.0 (20240811.2233)\n",
       " -->\n",
       "<!-- Pages: 1 -->\n",
       "<svg width=\"208pt\" height=\"290pt\"\n",
       " viewBox=\"0.00 0.00 208.34 289.50\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 285.5)\">\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-285.5 204.34,-285.5 204.34,4 -4,4\"/>\n",
       "<g id=\"clust1\" class=\"cluster\">\n",
       "<title>cluster_config_i</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"34.34,-156.5 34.34,-273.5 192.34,-273.5 192.34,-156.5 34.34,-156.5\"/>\n",
       "<text text-anchor=\"middle\" x=\"161.84\" y=\"-163.7\" font-family=\"Times,serif\" font-size=\"14.00\">config_i</text>\n",
       "</g>\n",
       "<g id=\"clust2\" class=\"cluster\">\n",
       "<title>cluster_param_i</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"42.34,-189 42.34,-265.5 112.34,-265.5 112.34,-189 42.34,-189\"/>\n",
       "<text text-anchor=\"middle\" x=\"81.84\" y=\"-196.2\" font-family=\"Times,serif\" font-size=\"14.00\">param_i</text>\n",
       "</g>\n",
       "<g id=\"clust3\" class=\"cluster\">\n",
       "<title>cluster_data</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"21.34,-8 21.34,-84.5 91.34,-84.5 91.34,-8 21.34,-8\"/>\n",
       "<text text-anchor=\"middle\" x=\"72.09\" y=\"-15.2\" font-family=\"Times,serif\" font-size=\"14.00\">data</text>\n",
       "</g>\n",
       "<!-- mu -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>mu</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" stroke-dasharray=\"5,2\" cx=\"103.34\" cy=\"-130.5\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"103.34\" y=\"-125.45\" font-family=\"Times,serif\" font-size=\"14.00\">mu</text>\n",
       "</g>\n",
       "<!-- obs -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>obs</title>\n",
       "<ellipse fill=\"grey\" stroke=\"black\" cx=\"56.34\" cy=\"-58.5\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"56.34\" y=\"-53.45\" font-family=\"Times,serif\" font-size=\"14.00\">obs</text>\n",
       "</g>\n",
       "<!-- mu&#45;&gt;obs -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>mu&#45;&gt;obs</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M92.68,-113.62C86.97,-105.11 79.8,-94.44 73.35,-84.82\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"76.29,-82.93 67.81,-76.58 70.48,-86.83 76.29,-82.93\"/>\n",
       "</g>\n",
       "<!-- scale -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>scale</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" cx=\"29.34\" cy=\"-130.5\" rx=\"29.34\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"29.34\" y=\"-125.45\" font-family=\"Times,serif\" font-size=\"14.00\">scale</text>\n",
       "</g>\n",
       "<!-- scale&#45;&gt;obs -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>scale&#45;&gt;obs</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M35.88,-112.55C38.86,-104.82 42.47,-95.46 45.83,-86.77\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"48.99,-88.29 49.32,-77.7 42.46,-85.77 48.99,-88.29\"/>\n",
       "</g>\n",
       "<!-- alpha -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>alpha</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" cx=\"153.34\" cy=\"-239.5\" rx=\"30.88\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"153.34\" y=\"-234.45\" font-family=\"Times,serif\" font-size=\"14.00\">alpha</text>\n",
       "</g>\n",
       "<!-- alpha&#45;&gt;mu -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>alpha&#45;&gt;mu</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M146.1,-221.73C139.65,-206.95 129.77,-184.66 117.03,-158.19\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"120.18,-156.67 112.67,-149.2 113.88,-159.73 120.18,-156.67\"/>\n",
       "</g>\n",
       "<!-- w -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>w</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" cx=\"77.34\" cy=\"-239.5\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"77.34\" y=\"-234.45\" font-family=\"Times,serif\" font-size=\"14.00\">w</text>\n",
       "</g>\n",
       "<!-- w&#45;&gt;mu -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>w&#45;&gt;mu</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M81.52,-221.31C85.56,-204.7 91.74,-179.24 96.5,-159.66\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"99.86,-160.66 98.82,-150.12 93.06,-159.01 99.86,-160.66\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x36ab792e0>"
      ]
     },
     "execution_count": 258,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numpyro.render_model(\n",
    "    model, \n",
    "    model_args=(\n",
    "        parameter_mask,\n",
    "        config_idx,\n",
    "        X,\n",
    "        y\n",
    "    ),\n",
    "    render_params=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 500/500 [00:03<00:00, 148.27it/s]\n"
     ]
    }
   ],
   "source": [
    "rng_key = random.PRNGKey(seed=42)\n",
    "rng_key, rng_subkey = random.split(key=rng_key)\n",
    "\n",
    "kernel = numpyro.infer.NUTS(model)\n",
    "\n",
    "mcmc = numpyro.infer.MCMC(\n",
    "    kernel, \n",
    "    num_warmup=200, \n",
    "    num_samples=300, \n",
    "    num_chains=4, \n",
    "    chain_method=\"vectorized\"\n",
    ")\n",
    "\n",
    "mcmc.run(\n",
    "    rng_subkey,\n",
    "    parameter_mask,\n",
    "    config_idx,\n",
    "    X,\n",
    "    y\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "  alpha[0]     -0.03      0.11     -0.04     -0.20      0.15    598.37      1.00\n",
      "  alpha[1]      0.05      0.12      0.05     -0.18      0.22    797.79      1.01\n",
      "     scale      0.52      0.02      0.52      0.50      0.55   1058.94      1.00\n",
      "    w[0,0]      2.51      0.01      2.51      2.49      2.53    837.47      1.00\n",
      "    w[0,1]      1.00      0.01      1.00      0.98      1.02   1350.56      1.00\n",
      "    w[1,0]      2.98      0.02      2.98      2.94      3.02    740.98      1.00\n",
      "    w[1,1]      0.49      0.03      0.49      0.45      0.53   1002.80      1.00\n",
      "    w[2,0]      0.11      5.30     -0.03     -7.42      9.98    893.63      1.00\n",
      "    w[2,1]      3.05      0.07      3.05      2.93      3.15    891.06      1.01\n",
      "\n",
      "Number of divergences: 0\n"
     ]
    }
   ],
   "source": [
    "# Parameter estimates for w[2, 0] are actually 0.0 once data is passed \n",
    "# through the program\n",
    "mcmc.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng_key, rng_subkey = random.split(key=rng_key)\n",
    "\n",
    "samples = mcmc.get_samples()\n",
    "\n",
    "predictive = Predictive(model, posterior_samples=samples)\n",
    "pps = predictive(\n",
    "    rng_subkey,\n",
    "    parameter_mask,\n",
    "    config_idx,\n",
    "    X,\n",
    "    None\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAE8CAYAAACIDWb/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJt0lEQVR4nO3deXxTVfo/8E/SJumSphvdW0qhULaCymYRQaCyuBWpsrkURB3ZFHCt31FAxwE3YBhZHEeBUZFFRRQVZK0DUmQZtGyFVoRCS1sgXZK2Sduc3x/9NZAm6UbS2+Xzfr3ykpxz782T02t5uPec58qEEAJERERE1KTkUgdARERE1BYxCSMiIiKSAJMwIiIiIgkwCSMiIiKSAJMwIiIiIgkwCSMiIiKSAJMwIiIiIgkwCSMiIiKSAJMwIiIiIgkwCSNqA9asWQOZTIY///xT6lCahEwmw/z5883vHf39//zzT8hkMqxZs8Yhx3OmvXv3QiaTYe/evU7/rPnz50Mmk1m0yWQyzJw50+mfDbS985xaPiZhRDZU/zKvfrm5uaFLly6YOXMmcnNzHf55JSUlmD9/fpP8Rels1X8RV788PDzQvXt3/PWvf0VRUZHU4TXIunXrsHTpUqnDMKtO/qpfCoUC7dq1w8CBA/Hqq6/iwoULDvusv//97/jmm28cdjxHas6xETWEjM+OJLK2Zs0aTJkyBW+88QaioqJQVlaGffv24dNPP0VkZCSOHz8ODw8Ph33elStXEBAQgHnz5llcwXGUyspKlJeXQ6VSWV2pcLT58+djwYIFWLlyJdRqNXQ6HX766Sds3rwZcXFx2L9/v9NjkMlkFmPZ2O9/33334fjx41ZXVoQQMBgMUCgUcHFxcWDktfvzzz8RFRWFiRMn4p577oHJZIJWq8WhQ4fw9ddfQyaT4eOPP8aECRPM+5hMJhiNRiiVSsjl9f93t1qtxkMPPdSgq30VFRWoqKiAm5ubuU0mk2HGjBn44IMP6n2cxsbWlOc5kSO4Sh0AUXM2evRo9O3bFwDw5JNPwt/fH4sXL8aWLVswceJEiaOrm16vh6enJ1xcXByaLJSUlNSZhD700ENo164dAOCZZ55BYmIivv76a6SmpiIuLq7Rx20MR3//6qujUrntttvw6KOPWrSdP38eI0aMQFJSErp164bevXsDAORyudNjrT7PXF1d4eoq3V8rjv45Ezkbb0cSNcCwYcMAAOfOnQNQ9S//N998E506dYJKpUKHDh3w6quvwmAwWOx3+PBhjBw5Eu3atYO7uzuioqLwxBNPAKi6uhEQEAAAWLBggflW041XxE6fPo2HHnoIfn5+cHNzQ9++ffHtt99afEb1LdSUlBRMnz4dgYGBCA8Pt+ireUVnxYoV6NGjB1QqFUJDQzFjxgwUFBRYbHPXXXehZ8+eOHLkCAYPHgwPDw+8+uqrNz12tR3XYDBg3rx5iI6OhkqlQkREBF566SWrcTUYDJgzZw4CAgLg5eWFBx54ABcvXrT6bHvf/8cff8SQIUPg5eUFjUaDfv36Yd26deb4vv/+e5w/f978M+nQoQMA6zlh7733HmQyGc6fP2/12cnJyVAqldBqtea2gwcPYtSoUfD29oaHhweGDBmC/fv3N3hMbxQZGYk1a9bAaDTinXfeMbfbmhN29uxZJCYmIjg4GG5ubggPD8eECRNQWFgIoCrJ1Ov1WLt2rfm7T548GcD1280nT57EpEmT4Ovri0GDBln02fL5558jJiYGbm5u6NOnD37++WeL/smTJ5vH90Y1j1lbbI44z0+ePImhQ4fCw8MDYWFhFmNJ5Gi8EkbUAJmZmQAAf39/AFVXx9auXYuHHnoIzz//PA4ePIiFCxfi1KlT2Lx5MwAgLy8PI0aMQEBAAF555RX4+Pjgzz//xNdffw0ACAgIwMqVKzFt2jQ8+OCDGDt2LACgV69eAIATJ07gjjvuQFhYGF555RV4enpi48aNGDNmDL766is8+OCDFjFOnz4dAQEBeP3116HX6+1+l+rbhvHx8Zg2bRrS09OxcuVKHDp0CPv374dCoTBve/XqVYwePRoTJkzAo48+iqCgoJseO3vHNZlMeOCBB7Bv3z48/fTT6NatG9LS0rBkyRKcOXPGYi7Qk08+ic8++wyTJk3CwIEDsXv3btx77731imfNmjV44okn0KNHDyQnJ8PHxwf/+9//sG3bNkyaNAn/93//h8LCQly8eBFLliwBUHUbzJZx48bhpZdewsaNG/Hiiy9a9G3cuBEjRoyAr68vAGD37t0YPXo0+vTpg3nz5kEul2P16tUYNmwY/vvf/6J///71HtOa4uLi0KlTJ+zYscPuNkajESNHjoTBYMCsWbMQHByMS5cuYevWrSgoKIC3tzc+/fRTPPnkk+jfvz+efvppAECnTp0sjvPwww+jc+fO+Pvf/466ZrWkpKRgw4YNePbZZ6FSqbBixQqMGjUKv/76K3r27Nmg71if2G7UkPNcq9Vi1KhRGDt2LMaNG4cvv/wSL7/8MmJjYzF69OgGxUlUL4KIrKxevVoAEDt37hT5+fkiKytLrF+/Xvj7+wt3d3dx8eJFcezYMQFAPPnkkxb7vvDCCwKA2L17txBCiM2bNwsA4tChQ3Y/Lz8/XwAQ8+bNs+obPny4iI2NFWVlZeY2k8kkBg4cKDp37mwV86BBg0RFRYXN73Pu3DkhhBB5eXlCqVSKESNGiMrKSvN2H3zwgQAgPvnkE3PbkCFDBACxatWqugdOCDFv3jwBQKSnp4v8/Hxx7tw58eGHHwqVSiWCgoKEXq+v9biffvqpkMvl4r///a9F+6pVqwQAsX//fiGEMI//9OnTLbabNGmS1VjW/P4FBQXCy8tLDBgwQJSWllrsbzKZzH++9957RWRkpNV3PHfunAAgVq9ebW6Li4sTffr0sdju119/FQDEf/7zH/OxO3fuLEaOHGnxOSUlJSIqKkrcfffdVp9l63Pfffddu9skJCQIAKKwsFAIIcSePXsEALFnzx4hhBD/+9//BACxadOmWj/L09NTJCUlWbVX/3wnTpxot+9GAAQAcfjwYXPb+fPnhZubm3jwwQfNbUlJSTbH2tYx7cXmiPO8+mclhBAGg0EEBweLxMREq88icgTejiSqRXx8PAICAhAREYEJEyZArVZj8+bNCAsLww8//AAAmDt3rsU+zz//PADg+++/BwD4+PgAALZu3Yry8vIGff61a9ewe/dujBs3DsXFxbhy5QquXLmCq1evYuTIkTh79iwuXbpksc9TTz1V57yYnTt3wmg0Yvbs2RaTtZ966iloNBpz7NVUKhWmTJnSoNhjYmIQEBCAqKgo/OUvf0F0dDS+//57izlfto67adMmdOvWDV27djV/3ytXrphvZ+7ZswcAzOP/7LPPWuw/e/bsOmPbsWMHiouL8corr1jNl2rshO7x48fjyJEj5it+ALBhwwaoVCokJCQAAI4dO4azZ89i0qRJuHr1qvm76fV6DB8+HD///DNMJlOjPr9a9dW64uJim/3e3t4AgO3bt6OkpKTRn/PMM8/Ue9u4uDj06dPH/L59+/ZISEjA9u3bUVlZ2egY6tLQ81ytVlvMtVMqlejfvz/++OMPp8VIbRtvRxLVYvny5ejSpQtcXV0RFBSEmJgY8y/z8+fPQy6XIzo62mKf4OBg+Pj4mOcHDRkyBImJiViwYAGWLFmCu+66C2PGjMGkSZOgUqlq/fyMjAwIIfDaa6/htddes7lNXl4ewsLCzO+joqLq/F7VscXExFi0K5VKdOzY0WpuU1hYGJRKZZ3HvdFXX30FjUYDhUKB8PBwm7eMbB337NmzOHXqlHmeXE15eXnm7yCXy62OW/M72VKdKDX0VlhtHn74YcydOxcbNmzAq6++CiEENm3ahNGjR0Oj0QCo+m4AkJSUZPc4hYWF5luXjaHT6QAAXl5eNvujoqIwd+5cLF68GJ9//jnuvPNOPPDAA3j00UfNCVp91Oc8q9a5c2erti5duqCkpAT5+fkIDg6u97EaoqHneXh4uFUS7uvri99//90p8RExCSOqRf/+/c2rI+2p68qJTCbDl19+idTUVHz33XfYvn07nnjiCbz//vtITU21O88IgPmqyAsvvICRI0fa3KZmEuju7l5rPI3RmGMOHjzYvDqyIcc1mUyIjY3F4sWLbe4TERHR4FiaQmhoKO68805s3LgRr776KlJTU3HhwgW8/fbb5m2qf57vvvsubrnlFpvHqe18qI/jx48jMDDQnPjZ8v7772Py5MnYsmULfvrpJzz77LNYuHAhUlNTzYs56uLo88ze/0fOvFJWk70ryIKVnMhJmIQRNVJkZCRMJhPOnj2Lbt26mdtzc3NRUFCAyMhIi+1vv/123H777Xjrrbewbt06PPLII1i/fj2efPJJu38BdezYEQCgUCgQHx/v0NgBID093fwZQNWk7XPnzjn0sxqqU6dO+O233zB8+PBaE9zq8c/MzLS40pGenl6vzwCqEpaaSeyNGnprcvz48Zg+fTrS09OxYcMGeHh44P7777f6XI1G45QxPnDgADIzM63KV9gSGxuL2NhY/PWvf8Uvv/yCO+64A6tWrcLf/vY3AI2/LWtL9RXAG505cwYeHh7mK56+vr5WKxYB2FxxWt/YmvN5TgSwRAVRo91zzz0AYFVRvfoKTvUqPa1Wa/Uv6eqrINUlF6rnSdX8SygwMBB33XUXPvzwQ+Tk5FjFkJ+f36jY4+PjoVQqsWzZMovYPv74YxQWFtZ7haEzjBs3DpcuXcJHH31k1VdaWmpe8Vm9Wm3ZsmUW29Snwv2IESPg5eWFhQsXoqyszKLvxvHw9PQ0l22oj8TERLi4uOCLL77Apk2bcN9998HT09Pc36dPH3Tq1Anvvfee+bbhjRr78wSqkpXJkydDqVRardC8UVFRESoqKizaYmNjIZfLLUqAeHp62kyKGuPAgQM4evSo+X1WVha2bNmCESNGmK8+derUCYWFhRa3/nJycsyrjG9U39ia83lOBPBKGFGj9e7dG0lJSfjXv/6FgoICDBkyBL/++ivWrl2LMWPGYOjQoQCAtWvXYsWKFXjwwQfRqVMnFBcX46OPPoJGozEncu7u7ujevTs2bNiALl26wM/PDz179kTPnj2xfPlyDBo0CLGxsXjqqafQsWNH5Obm4sCBA7h48SJ+++23BsceEBCA5ORkLFiwAKNGjcIDDzyA9PR0rFixAv369avXlRRneeyxx7Bx40Y888wz2LNnD+644w5UVlbi9OnT2LhxI7Zv346+ffvilltuwcSJE7FixQoUFhZi4MCB2LVrFzIyMur8DI1GgyVLluDJJ59Ev379zPWufvvtN5SUlGDt2rUAqpKmDRs2YO7cuejXrx/UarXFla2aAgMDMXToUCxevBjFxcUYP368Rb9cLse///1vjB49Gj169MCUKVMQFhaGS5cuYc+ePdBoNPjuu+/qjP/o0aP47LPPYDKZUFBQgEOHDuGrr76CTCbDp59+ai5vYsvu3bsxc+ZMPPzww+jSpQsqKirw6aefwsXFBYmJiebt+vTpg507d2Lx4sUIDQ1FVFQUBgwYUGdstvTs2RMjR460KFEBVNXFqzZhwgS8/PLLePDBB/Hss8+ipKQEK1euRJcuXSwSuIbE1pzPcyIALFFBZEv1UvfaykoIIUR5eblYsGCBiIqKEgqFQkRERIjk5GSLchJHjx4VEydOFO3btxcqlUoEBgaK++67z2LJvhBC/PLLL6JPnz5CqVRalVjIzMwUjz/+uAgODhYKhUKEhYWJ++67T3z55Zf1irnm0v1qH3zwgejatatQKBQiKChITJs2TWi1WotthgwZInr06FHHiF1XXVIgPz+/1u1qO67RaBRvv/226NGjh1CpVMLX11f06dNHLFiwwFx6QQghSktLxbPPPiv8/f2Fp6enuP/++0VWVladJSqqffvtt2LgwIHC3d1daDQa0b9/f/HFF1+Y+3U6nZg0aZLw8fERAMwlFGyVqKj20UcfCQDCy8vLqvxFtf/9739i7Nixwt/fX6hUKhEZGSnGjRsndu3aVeuYVX9u9cvV1VX4+fmJAQMGiOTkZHH+/HmrfWqWqPjjjz/EE088ITp16iTc3NyEn5+fGDp0qNi5c6fFfqdPnxaDBw8W7u7uAoC5JERtP197JSpmzJghPvvsM9G5c2ehUqnErbfeao7nRj/99JPo2bOnUCqVIiYmRnz22Wc2j2kvNmec5/ZKZxA5Ap8dSURERCQBzgkjIiIikgCTMCIiIiIJMAkjIiIikgCTMCIiIiIJMAkjIiIikgCTMCIiIiIJtPpirSaTCdnZ2fDy8nLoYziIiIiIahJCoLi4GKGhoZDLa7/W1eqTsOzs7Gb7wF8iIiJqnbKyshAeHl7rNq0+CfPy8gJQNRgajUbiaIiIiKg1KyoqQkREhDn/qE2rT8Kqb0FqNBomYURERNQk6jMFihPziYiIiCTAJIyIiIhIAkzCiIiIiCTAJIyIiIhIAkzCiIiIiCTAJIyIiIhIAkzCiIiIiCTQ6uuEUduUk5MDrVbboH18fX0REhLipIiIiIgsMQmjVicnJwcxXTqjWKdv0H5eak+knznLRIyIiJoEkzBqdbRaLYp1emx5ti+iAz3qtU9GXgkSlh2GVqtlEkZERE2CSRi1WtGBHugeVvezu4iIiKTAiflEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEREREEuDqSGpTDAYDKioqrNpLSksBABkZGRBCWPX7+voiNDTU6fEREVHbwSSM2gyDwYCDBw+isrLSqu/PAhMAICEhwea+nmovnEk/zUSMiIgchkkYtRkVFRWorKyEd1hHuChVFn3a/FIAJzFo+ttQB4RZ9OnyLmLfyleg1WqZhBERkcMwCaM2x0WpgqvS3aLNVVF1C1IdEA7v0CgpwiIiojaGE/OJiIiIJMAkjIiIiEgCTMKIiIiIJMAkjIiIiEgCTMKIiIiIJMAkjIiIiEgCTMKIiIiIJMAkjIiIiEgCTMKIiIiIJMAkjIiIiEgCTMKIiIiIJMAkjIiIiEgCTMKIiIiIJNBskrBFixZBJpNh9uzZ5raysjLMmDED/v7+UKvVSExMRG5urnRBEhERETlIs0jCDh06hA8//BC9evWyaJ8zZw6+++47bNq0CSkpKcjOzsbYsWMlipKIiIjIcSRPwnQ6HR555BF89NFH8PX1NbcXFhbi448/xuLFizFs2DD06dMHq1evxi+//ILU1FS7xzMYDCgqKrJ4ERERETU3kidhM2bMwL333ov4+HiL9iNHjqC8vNyivWvXrmjfvj0OHDhg93gLFy6Et7e3+RUREeG02ImIiIgaS9IkbP369Th69CgWLlxo1Xf58mUolUr4+PhYtAcFBeHy5ct2j5mcnIzCwkLzKysry9FhExEREd00V6k+OCsrC8899xx27NgBNzc3hx1XpVJBpVI57HhEREREziDZlbAjR44gLy8Pt912G1xdXeHq6oqUlBQsW7YMrq6uCAoKgtFoREFBgcV+ubm5CA4OliZoIiIiIgeR7ErY8OHDkZaWZtE2ZcoUdO3aFS+//DIiIiKgUCiwa9cuJCYmAgDS09Nx4cIFxMXFSREyERERkcNIloR5eXmhZ8+eFm2enp7w9/c3t0+dOhVz586Fn58fNBoNZs2ahbi4ONx+++1ShExERETkMJIlYfWxZMkSyOVyJCYmwmAwYOTIkVixYoXUYRERERHdtGaVhO3du9fivZubG5YvX47ly5dLExARERGRk0heJ4yIiIioLWpWV8KIpKbLv2ij7RIAIDMzEzKZzKLP19cXISEhTRIbERG1LkzCiABc05dDLgP2rXjZ7jYJCQlWbV5qT6SfOctEjIiIGoxJGBEAXVkFTAJYOiEaUUEai76K8jIUXvwDsbGx8HB3N7dn5JUgYdlhaLVaJmFERNRgTMKIbhDhp0KnQA+LtgqjDNd0cnQL8YSnp6dEkRERUWvDiflEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEREREEmASRkRERCQBJmFEREREEnCVOgCixsrOzoZWq7Vqz8zMBACUlJZCr7/+74yS0tKb+rya+1e/3717NzIyMqy2NxqNUCqVVu0ajQYBAQE2P8PX1xehoaE3FScREbUMTMKoRcrOzkaXmK7Q64rtbpOWlobiLOuLvcIkGvRZpooKQAYcT0uzaD92uRJyGTBr1qwGHU8uA+yF4Kn2wpn000zEiIjaACZh1CJptVrodcUYNG0R1IHhFn26/EvYt+JleId3hF+Au7ndqNdBl5cFIRqWhAlTJSAAr5AoKNzczO1yYwFMIhOLH4pEZKCnxT7lJSUouXYZngFhcFWpzO1Z1wyYvT4Tg6a/DXVAmGXceRexb+Ur0Gq1TMKIiNoAJmHUoqkDw+Ed2rFGqwwA4Kpwg6vyehJWaTTc1Ge5KJQWx3NxLQEARAaqERPma7GtQSdHoUkOnxANlO7XEzRXRdU+6oBweIdG3VQ8RETUsnFiPhEREZEEmIQRERERSYBJGBEREZEEmIQRERERSYBJGBEREZEEmIQRERERSYBJGBEREZEEmIQRERERSYBJGBEREZEEmIQRERERSYBJGBEREZEEmIQRERERSUDSJGzlypXo1asXNBoNNBoN4uLi8OOPP5r7y8rKMGPGDPj7+0OtViMxMRG5ubkSRkxERETkGJImYeHh4Vi0aBGOHDmCw4cPY9iwYUhISMCJEycAAHPmzMF3332HTZs2ISUlBdnZ2Rg7dqyUIRMRERE5hKuUH37//fdbvH/rrbewcuVKpKamIjw8HB9//DHWrVuHYcOGAQBWr16Nbt26ITU1FbfffrsUIRMRERE5hKRJ2I0qKyuxadMm6PV6xMXF4ciRIygvL0d8fLx5m65du6J9+/Y4cOCA3STMYDDAYDCY3xcVFTk9dqKG0uVftNF2CQCQmZkJmUxm1e/r64uQkBCnx0ZERE1D8iQsLS0NcXFxKCsrg1qtxubNm9G9e3ccO3YMSqUSPj4+FtsHBQXh8uXLdo+3cOFCLFiwwMlREzXONX055DJg34qX7W6TkJBgs91L7Yn0M2eZiBERtRKSJ2ExMTE4duwYCgsL8eWXXyIpKQkpKSmNPl5ycjLmzp1rfl9UVISIiAhHhEp003RlFTAJYOmEaEQFaSz6KsrLUHjxD8TGxsLD3d2iLyOvBAnLDkOr1TIJIyJqJSRPwpRKJaKjowEAffr0waFDh/CPf/wD48ePh9FoREFBgcXVsNzcXAQHB9s9nkqlgkqlcnbYRDclwk+FToEeFm0VRhmu6eToFuIJT09PiSIjIqKm0uzqhJlMJhgMBvTp0wcKhQK7du0y96Wnp+PChQuIi4uTMEIiIiKimyfplbDk5GSMHj0a7du3R3FxMdatW4e9e/di+/bt8Pb2xtSpUzF37lz4+flBo9Fg1qxZiIuL48pIIiIiavEkTcLy8vLw+OOPIycnB97e3ujVqxe2b9+Ou+++GwCwZMkSyOVyJCYmwmAwYOTIkVixYoWUIRMRERE5hKRJ2Mcff1xrv5ubG5YvX47ly5c3UURERERETaPZzQkjIiIiaguYhBERERFJgEkYERERkQQalYR17NgRV69etWovKChAx44dbzooIiIiotauUUnYn3/+icrKSqt2g8GAS5cu3XRQRERERK1dg1ZHfvvtt+Y/V9fyqlZZWYldu3ahQ4cODguOiIiIqLVqUBI2ZswYAIBMJkNSUpJFn0KhQIcOHfD+++87LDgiIiKi1qpBSZjJZAIAREVF4dChQ2jXrp1TgiIiIiJq7RpVrPXcuXOOjoOIiIioTWl0xfxdu3Zh165dyMvLM18hq/bJJ5/cdGBERERErVmjkrAFCxbgjTfeQN++fRESEgKZTObouIiIiIhatUYlYatWrcKaNWvw2GOPOToeIiIiojahUXXCjEYjBg4c6OhYiIiIiNqMRiVhTz75JNatW+foWIiIiIjajEbdjiwrK8O//vUv7Ny5E7169YJCobDoX7x4sUOCIyIiImqtGpWE/f7777jlllsAAMePH7fo4yR9IiIioro1Kgnbs2ePo+OgViI7OxtarbbB+/n6+iI0NNQJERERETVPja4TRlRTdnY2usR0hV5X3OB9PdVeOJN+mokYERG1GY1KwoYOHVrrbcfdu3c3OiBqubRaLfS6YgyatgjqwPB676fLu4h9K1+BVqtlEkZERG1Go5Kw6vlg1crLy3Hs2DEcP37c6sHe1PaoA8PhHdpR6jCIiIiatUYlYUuWLLHZPn/+fOh0upsKiIiIiKgtcOicsEcffRT9+/fHe++958jDUhuWk5Njc6J/ZmYmAECXfwmA5a1xXf7FpgiNiIjopjg0CTtw4ADc3NwceUhqw3JychDTpTOKdXq72+xb8bLdPlNFuTPCIiIicohGJWFjx461eC+EQE5ODg4fPozXXnvNIYERabVaFOv02PJsX0QHelj0lZSWIi0tDd7hHeGqsEz8D2YW4O9b/4BJmJoyXCIiogZpVBLm7e1t8V4ulyMmJgZvvPEGRowY4ZDAiKpFB3qge5iXRZteL0dxlhx+Ae5wVbpb9F24WtqU4RERETVKo5Kw1atXOzoOIiIiojblpuaEHTlyBKdOnQIA9OjRA7feeqtDgiJqy0pKra/kVbdlZGRACGHVzycOEBG1PI1KwvLy8jBhwgTs3bsXPj4+AICCggIMHToU69evR0BAgCNjJGoTTBUVgAw4npZm1fdnQdX8toSEBJv78okDREQtT6OSsFmzZqG4uBgnTpxAt27dAAAnT55EUlISnn32WXzxxRcODZKoLRCmSkAAXiFRUNRYZazNLwVwEoOmvw11QJhFH584QETUMjUqCdu2bRt27txpTsAAoHv37li+fDkn5hPdJBeF0mqxgaui6hakOiAc3qFRUoRFREQOJm/MTiaTCQqFwqpdoVDAZGJZACIiIqK6NOpK2LBhw/Dcc8/hiy++MN/+uHTpEubMmYPhw4c7NEBqOzIyMizeV1fFLykthV5v+e8FW5PXiYiIWpJGJWEffPABHnjgAXTo0AEREREAgKysLPTs2ROfffaZQwOk1q+sWAvIZBgzZozN/rS0NBRn2b5oK0zWKwWJiIhagkYlYRERETh69Ch27tyJ06dPAwC6deuG+Ph4hwZHbUN5qR4QAn2fWIB2EZ3M7br8S9i34mV4h3eEX4DlHCmjXgddXpbNcg1EREQtQYOSsN27d2PmzJlITU2FRqPB3XffjbvvvhsAUFhYiB49emDVqlW48847nRIstW6e7cLgHdrxhpaqB3O7KtysJqpXGg1NGBkREZHjNSgJW7p0KZ566iloNBqrPm9vb/zlL3/B4sWLmYSRTaWFV2DUF1u1l1zLBQDor2Sj0P16sqXLv9hksRERETW1BiVhv/32G95++227/SNGjMB7771300FR61NaeAXbXh8HY1mZ3W0Of/K6zXZTRbmzwiIiIpJMg5Kw3Nxcm6UpzAdzdUV+fn69j7dw4UJ8/fXXOH36NNzd3TFw4EC8/fbbiImJMW9TVlaG559/HuvXr4fBYMDIkSOxYsUKBAUFNSR0kphRXwxjWRn++Wg3tPevMb+rpAjFuRehCekAhZuHuf1gZgH+vvUPmATLnhARUevToDphYWFhOH78uN3+33//HSEhIfU+XkpKCmbMmIHU1FTs2LED5eXlGDFiBPR6vXmbOXPm4LvvvsOmTZuQkpKC7OxsjB07tiFhUzPS3t8dnQI9LF4d27mhg48cHQMs+0J8VFKHS0RE5DQNuhJ2zz334LXXXsOoUaPgVuOxKqWlpZg3bx7uu+++eh9v27ZtFu/XrFmDwMBAHDlyBIMHD0ZhYSE+/vhjrFu3DsOGDQMArF69Gt26dUNqaipuv/32hoRPRERE1Gw0KAn761//iq+//hpdunTBzJkzzbcNT58+jeXLl6OyshL/93//1+hgCgsLAQB+fn4AgCNHjqC8vNyi9EXXrl3Rvn17HDhwwGYSZjAYYDBcXzlXVFTU6HjIvpycHGi1Wou26uKquvxLqF7ZWI2T7J2vZrHbuvj6+vJZk0REEmpQEhYUFIRffvkF06ZNQ3JysrlGk0wmw8iRI7F8+fJGz9UymUyYPXs27rjjDvTs2RMAcPnyZSiVSvj4+FjFcfnyZZvHWbhwIRYsWNCoGKh+cnJyENOlM4p1epv9+1a8bHdfTrJ3vLqK3drjqfbCmfTTTMSIiCTS4GKtkZGR+OGHH6DVapGRkQEhBDp37gxfX9+bCmTGjBk4fvw49u3bd1PHSU5Oxty5c83vi4qKzFX9yTG0Wi2KdXpsebYvogOvT6QvKS1FWloavMM7wlVhebuak+ydx16x29ro8i5i38pXoNVqmYQREUmkURXzgapbGf369XNIEDNnzsTWrVvx888/Izw83NweHBwMo9GIgoICi6thubm5CA4OtnkslUoFlYoTuptCdKAHuod5md/r9XIUZ8nhF+BuVVz1wlU+69HZrIvdEhFRc9ag1ZGOJoTAzJkzsXnzZuzevRtRUVEW/X369IFCocCuXbvMbenp6bhw4QLi4uKaOlwiIiIih2n0lTBHmDFjBtatW4ctW7bAy8vLPM/L29sb7u7u8Pb2xtSpUzF37lz4+flBo9Fg1qxZiIuL48pIapNsLXCw98QBADBVlkPuYl3br2rxRNViCpnMchGFr69vg0rNEBFR40iahK1cuRIAcNddd1m0r169GpMnTwYALFmyBHK5HImJiRbFWonakmv6cshltS96sPXEAbkMMNXyjPOEhASrNi+1J9LPnGUiRkTkZJImYdWrK2vj5uaG5cuXY/ny5U0QEVHzpCurgEkASydEIyrI8tmtdT1xwNY+FeVlKLz4B2JjY+Fxw9WzjLwSJCw7DK1WyySMiMjJJE3CiKhhIvxU6HTDilQAMOiMKDTI4RPgDqX79b7qxRC29qkwynBNJ0e3EE94eno6P3AiIrIi6cR8IiIioraKV8LaOFuV7+vS0MrsREREZI1JWBtWV+X7uhjLjQ6OiIiIqO1gEtaG2at8X5fdp65i1ucnUFFR6cToiIiIWjcmYWRV+b4uGXmNu3JGRERE13FiPhEREZEEmIQRERERSYBJGBEREZEEmIQRERERSYBJGBEREZEEuDqyjcjOzrYqypqZmQkAKCkthV5vnY+7urpCpVI1SXwkjZLSUpvvMzIybD7b1dfXF6GhoU0SW1vSmKLJQNXPg8/4JGq5mIS1AdnZ2egS0xV6XbHN/rS0NBRnWSdhLi4uGDBgABOxVshUUQHIgONpaRbtfxaYAAAJCQk29/NUe+FM+mkmYg50M0WTvdSeSD9zlokYUQvFJKwN0Gq10OuKMWjaIqgDw83tuvxL2LfiZXiHd4RfgLvFPpVGAwov/YGKigomYa2QMFUCAvAKiYLCzc3crs0vBXASg6a/DXVAmMU+uryL2LfyFWi1WiZhDtTYoskZeSVIWHYYWq2WSRhRC8UkrA1RB4bDO7TjDS0yAICrwg2uSnfbO1Gr5qJQWvzsXRVVtyDVAeHwDo2SKqw2qaFFk4mo5ePEfCIiIiIJ8EoY1armxG0AKCszAAAMhjLo9fpatyUCbC8MqYu9RQC1TWLPy8tDUVGRVXt5eTkUCoXNfTQaDQIDA+3GYO9WX2O+k8FgsLq939QLZBqzCIALAIicg0kY2WRv4jYAnMmqqPrvmbMw5mVa9QuT9ao6arvqWhhij61FAI2dxC6XAY05Le1NfG/sd4JMDgiTza6mWCDT2PHjAgAi52ASRjbZm7gNAF7GAgCZ8AqKgF+oxtxu1Ougy8uyWdqA2i57C0NqY28RQG2T2EtKS5GWlgavoHDIFUpz+6FzxVj0YxYWPxSJyEBPi31M5UYU515EbGwsPNwt50XWNvG9Md8p9/QR/G/DEvR9YgHaRXS6/l2bcIFMYxYBcAEAkfMwCaNa1Zy4DQAuriUAAHmNvkqjoUljo5bFemFI49maxK7Xy1GcJYdfqLfFeZldVPWPgshANWLCfC32qTCW4pohG91CPOHpaZmg1UdDvlNx3kUAgGe7MMkXyHARAFHzwIn5RERERBJgEkZEREQkASZhRERERBJgEkZEREQkASZhRERERBJgEkZEREQkAZaoaGFqq9Kdn59vs1r4hQsXAFTVI6peDl/1/qJTYiRylIyMDIv3tVWXb81PbLD13arbMjIybNbms/fEASJqPpiEtSB1Vemuqyr4vhUv22w3VZQ7Ijwihykr1gIyGcaMGWOz3151eaB1PbGhtidX/FlQVXk/ISHB5r62njhARM0Lk7AWpLYq3dVVt5dO6IQIP8vK2uUlJSi5dhmakA5QuF2vkn0wswB/3/oHTHYeo0IklfJSPSBEg6rLt8YnNtT25AptfimAkxg0/W2oA8Is+uw9cYCImhcmYS2Q7SrdVbcZo4K80anG40gMOjkKTXL4BLhD6X6978LV1nv7hlqHhlSXb81PbLD15ApXRVWyqQ4Ih3dolBRhEdFNYhJGRFZszResmlNYNS9LJpNZ9RuNRiiVSqv26nlcNeckAoCpshxyF4XVPiXXcgEA+ivZKLzhmY6cx9gwDZlTBwCurq4OeUYl1V9t83xrwzl/rQOTMCIyu6Yvh1xmf/4gYH8OkosMqGzgnMS65jEe/uR1m+2cx1i7xs6pc3FxwYABA5iINZG65vnWhnP+WgcmYURkpiurgEkASydEIypIY9FXUV6Gwot/IDY2Fh7ulrfGdp+6ilmfn8CGp3uiZ3s/i76S0lKkpaXBO7wjXBXX5zVVz0m09VnGkiIU517kPMZGasycukqjAYWX/kBFRQWTsCZS2zzf2nDOX+vBJIyIrET4qazmFlYYZbimk6NbiCc8PT0t+jLy9ACAjgHu6B7mZdGn18tRnCWHX4C7xbym6jmJtj7LoDOi0MB5jDerIXPqSDq25/lSW8BirUREREQS4JWwZionJ8dqsmZtE5w5YZnIOWpObgda92KDmoVhm7oorK3ffXXx9fVFSEiIw2IgaipMwpqhnJwcxHTpjGKd3mZ/bZOmOWGZyDHyiw2Qy+wvRABa12IDe4Vhm7IobF2/++zxUnsi/cxZJmLU4kiahP3888949913ceTIEeTk5GDz5s0Wq3mEEJg3bx4++ugjFBQU4I477sDKlSvRuXNn6YJuAlqtFsU6PbY82xfRN8yVsTfBGeCEZSJHKyytWqTQVhYb2CsM25RFYe397qtNRl4JEpYdhlarZRJGLY6kSZher0fv3r3xxBNPYOzYsVb977zzDpYtW4a1a9ciKioKr732GkaOHImTJ0/CrUb16NYoOtDDYpKzvQnOACcsEzlLW1tsULMwrBRFYWv+7iNqrSRNwkaPHo3Ro0fb7BNCYOnSpfjrX/9qvgz+n//8B0FBQfjmm28wYcKEpgyViIiIyKGa7erIc+fO4fLly4iPjze3eXt7Y8CAAThw4IDd/QwGA4qKiixeRERERM1Ns03CLl++DAAICgqyaA8KCjL32bJw4UJ4e3ubXxEREU6Nk4iIiKgxmm0S1ljJyckoLCw0v7KysqQOiYiIiMhKs03CgoODAQC5ubkW7bm5ueY+W1QqFTQajcWLiIiIqLlptklYVFQUgoODsWvXLnNbUVERDh48iLi4OAkjIyIiIrp5kq6O1Ol0FtWoz507h2PHjsHPzw/t27fH7Nmz8be//Q2dO3c2l6gIDQ21qCVGRE2rZkV1ACgrMwAADIYy6PX6Orcnx7FVcb+pq/PbeqoAAOTn59tdHKXRaBAYGFiv47QU2dnZNqv92xuHCxcuALD95AWlpxfcvds5JU5qPiRNwg4fPoyhQ4ea38+dOxcAkJSUhDVr1uCll16CXq/H008/jYKCAgwaNAjbtm1rEzXCiJobexXVAeBMVkXVf8+chTEv0+b+orYy8tRg1/TlkMtqf4KGs6vzlxVrAZnM7j+M63p6gD3GcuPNBSaB7OxsdInpCr2u2KqvrnGw9TNUurlh1BsbmYi1cpImYXfddZfNZ5FVk8lkeOONN/DGG280YVREZIu9iuoA4GUsAJAJr6AI+IXWqBSv10GXl1Xr/+vUcLqyqor+UlbnLy/VA0Kg7xML0C6ik2V8+Zewb8XLWDqhEyL8VBZ9pnIjinMvIjY2Fh43XKnbfeoqZn1+AhUVlQ6JrylptVrodcUYNG0R1IHh5vbaxqG8pAQl1y5b/ZwuXC3FrM9OwagvZhLWyvHZkUTUIDUrqgOAi2sJAEBuo6/SaGiy2Nqi5lCd37NdGLxDO9Zorbq9FhXkbRVfhbEU1wzZ6BbiCU9PT3N7Rl7DnhnZHKkDw2uMhf1xMOjkKDRZ/5yo7WASRkRtHue5tQwGgwEVFRUWbdU/i4yMDJtXW319fR3yXEsp2JvDVzWHDMjMzIRMZjmXzNfXl8/QbEGYhBFRm8V5bi2HwWDAwYMHUVlpeavyz4KqW6vVj7eryVPthTPpp1tUIlaf+X6A7e/spfZE+pmzTMRaCCZhRNRmcZ5by1FRUYHKykp4h3WEi/L63CptfimAkxg0/W2oA8Is9tHlXcS+la9Aq9W2qCSstvl+AFBRXobCi39YzanLyCtBwrLD0Gq1TMJaCCZhRNTmcZ5by+GiVFn8PFwVVYmwOiAc3qFRUoXlFLbm+wFAhVGGazq51Zw6anmabbFWIiIiotaMV8IkZK+wX2Zm1fyTktJS6PXX82ROBiYiqTSkKKy97a32r/E7rSkXQ+Tk5Nj8/VstLy/PZoFVjUaDgIAAq/aWXmiWpMEkTCK1FfarlpaWhuIs64uVnAxMRE3lZorCArYLw9pbENFUiyFycnIQ06UzinUNL4lRV+HVcqNjCuFS28AkTCL2CvsB14v7eYd3hF/A9X9dcjIwETW1xhSFBWovDGtvQURTLYbQarUo1umx5dm+iLYx56qktBRpaWnwCgqHXKE0t2ddM2D2+kybiwByTx/B/zYsQUVlyys0S9JhEiYx68J+QHVxP1eFm8UEVE4GJiKpNKQoLFC/wrA1F0Q09WKI6EAPdA/zsmrX6+UozpLDL9S7xiKAqvhsLQIoznPOczmpdePEfCIiIiIJ8EoYERG1eLYWAtirLH+zk+gbskihPgsUqO1iEkZERC1WfRYO2Kumbyw3Ovyz7C1SsLVAgYhJGBERtVi1LRywV1l+96mrmPX5CVRUNGwSfWMWKdS2QIGISRgREbV4thYO2Kssn5HX8NIUdX2WvUUK9VmgQG0XJ+YTERERSYBXwpzMXlXm6qr4VRNHZRZ9nMhJROQ4DanMb2t7ImdhEuZE9anKXNsET07kJCJqvJupzA/w6STkfEzCnKi2qszVFZm9wzvCVeFm0ceJnEREN68xlfkBPp2Emg6TsCZgqyqzuSJzgLtVZWhO5CQicpyGVOYH+HQSajqcmE9EREQkAV4JIyIiagOys7NtLhSrja+vL0JDQ50UETEJIyIiauWys7PRJaYr9LriBu3nqfbCmfTTTMSchEkYERFRK6fVaqHXFWPQtEVQB4bXax9d3kXsW/kKtFotkzAnYRJGRETURqgDw+Ed2lHqMOj/YxJGRETUhpUWXoFRb32bsqqYeFVxcZnMsqi40WiEUqm0e8y8vDwUFRVZtGk0GgQEBNjdpy3OP2MSRkRE1EaVFl7BttfHwVhWZnebhIQEqzYXGVDZwDJqchlQW/3btjj/jEkYERFRG2XUF8NYVoZ/PtoN7f0ta6ZVlJeh8OIfiI2NhYf79b7dp65i1ucnsOHpnujZ3s/qmNXFyL2CwiFXVF0ty7pmwOz1mRg0/W2oA8Ks9mmr88+YhBEREbVx7f3d0anGk10qjDJc08nRLcQTnp6e5vaMvKpH8XUMcLcqRA7cUIw81NtcDNdVUVUgVx0QDu/QKGd9jRaHxVqJiIiIJMArYQ5iqwheZmbVg2FLSkuh11vmuyWlfDQRERE1fzX/viorq3qsk8FQBr1eX+f2N9LlX7TT3vhFALb4+voiJCSkQftIgUmYA9RVBC8tLQ3FWbYvOoraZikSERFJxFRRAciA42lpFu1nsiqq/nvmLIx5mXb3v/Hvt2v6cshlwL4VL9f6mY5aBOCl9kT6mbPNPhFjEuYA9org6fIvYd+Kl+Ed3hF+AZYTHo16HXR5WRCCSRgRETU/wlQJCMArJAoKNzdzu5exAEAmvIIi4BeqsdrP1t9vurIKmASwdEI0ooKs92nsIgBbMvJKkLDsMLRaLZOwtsS6CF7VJVVXhZt5cmK1SqOhCSMjIiJqHBeF0uLvMBfXqkn28hrt1Wr7+y3CT2W1AABo/CKAlo4T84mIiIgkwCthRERELVDNCfDV73fv3o2MjAyLvgsXLgCongB/feK7vYnyUmnoIgCTyQS53PbCt4yMDLtTfppLdX4mYURERC2IvQnzxy5XQi4DZs2aZXdfexPjTRXlDo2xoRq9CEAGoEae9WeBCYDtSf7Vmkt1/haRhC1fvhzvvvsuLl++jN69e+Of//wn+vfvL3VYRERETc7ehHm5sQAmkYnFD0UiMtDTYp/ykhKUXLsMTUgHKNyuz8k6mFmAv2/9AyZharL4bWnMIoDqBQA199HmlwI42SKq8zf7JGzDhg2YO3cuVq1ahQEDBmDp0qUYOXIk0tPTERgYKHV4REREkrA3YT4yUI2YMF+LbQ06OQpNcvgEuEPpfj0Ju3C1edWsbMgigOoFADX3cVVUXRprCdX5m/3E/MWLF+Opp57ClClT0L17d6xatQoeHh745JNPpA6NiIiIqNGa9ZUwo9GII0eOIDk52dwml8sRHx+PAwcO2NzHYDDAYLi+PLawsBAAUFRU5LQ4dTodAKDgYgYqDCXX26/kAABOXdSiqMaEwvKSUpRqK+EhCuCqKrHo++Ny1fHSLxWhxFhxU/vUtl9j9mF8Nxefo78T42ub34nxOW8fxuec+JryO13SGgEABZcyUGGwnsyvy8+u+q9O55TcoPqY9aoDKpqxS5cuCQDil19+sWh/8cUXRf/+/W3uM2/ePIGqaXp88cUXX3zxxRdfkryysrLqzHOa9ZWwxkhOTsbcuXPN700mE65duwZ/f3+r51G1VkVFRYiIiEBWVhY0GuvKxG0Jx+I6jsV1HIvrOBbXcSyqcByua8xYCCFQXFxcr0n/zToJa9euHVxcXJCbm2vRnpubi+DgYJv7qFQqqFQqizYfHx9nhdisaTSaNv8/UDWOxXUci+s4FtdxLK7jWFThOFzX0LHw9vau13bNemK+UqlEnz59sGvXLnObyWTCrl27EBcXJ2FkRERERDenWV8JA4C5c+ciKSkJffv2Rf/+/bF06VLo9XpMmTJF6tCIiIiIGq3ZJ2Hjx49Hfn4+Xn/9dVy+fBm33HILtm3bhqCgIKlDa7ZUKhXmzZtndVu2LeJYXMexuI5jcR3H4jqORRWOw3XOHguZEPVZQ0lEREREjtSs54QRERERtVZMwoiIiIgkwCSMiIiISAJMwoiIiIgkwCSsFZk/fz5kMpnFq2vXrlKH1SR+/vln3H///QgNDYVMJsM333xj0S+EwOuvv46QkBC4u7sjPj4eZ8+elSZYJ6trLCZPnmx1nowaNUqaYJ1o4cKF6NevH7y8vBAYGIgxY8YgPT3dYpuysjLMmDED/v7+UKvVSExMtCoO3RrUZyzuuusuq/PimWeekShi51m5ciV69eplLr4ZFxeHH3/80dzfVs4JoO6xaCvnRE2LFi2CTCbD7NmzzW3OOi+YhLUyPXr0QE5Ojvm1b98+qUNqEnq9Hr1798by5ctt9r/zzjtYtmwZVq1ahYMHD8LT0xMjR45EWVlZE0fqfHWNBQCMGjXK4jz54osvmjDCppGSkoIZM2YgNTUVO3bsQHl5OUaMGAG9/voDfefMmYPvvvsOmzZtQkpKCrKzszF27FgJo3aO+owFADz11FMW58U777wjUcTOEx4ejkWLFuHIkSM4fPgwhg0bhoSEBJw4cQJA2zkngLrHAmgb58SNDh06hA8//BC9evWyaHfaeXHTT9mmZmPevHmid+/eUochOQBi8+bN5vcmk0kEBweLd99919xWUFAgVCqV+OKLLySIsOnUHAshhEhKShIJCQmSxCOlvLw8AUCkpKQIIarOAYVCITZt2mTe5tSpUwKAOHDggFRhNomaYyGEEEOGDBHPPfecdEFJyNfXV/z73/9u0+dEteqxEKLtnRPFxcWic+fOYseOHRbf3ZnnBa+EtTJnz55FaGgoOnbsiEceeQQXLlyQOiTJnTt3DpcvX0Z8fLy5zdvbGwMGDMCBAwckjEw6e/fuRWBgIGJiYjBt2jRcvXpV6pCcrrCwEADg5+cHADhy5AjKy8stzouuXbuiffv2rf68qDkW1T7//HO0a9cOPXv2RHJyMkpKSqQIr8lUVlZi/fr10Ov1iIuLa9PnRM2xqNaWzokZM2bg3nvvtfj5A879XdHsK+ZT/Q0YMABr1qxBTEwMcnJysGDBAtx55504fvw4vLy8pA5PMpcvXwYAq6csBAUFmfvaklGjRmHs2LGIiopCZmYmXn31VYwePRoHDhyAi4uL1OE5hclkwuzZs3HHHXegZ8+eAKrOC6VSCR8fH4ttW/t5YWssAGDSpEmIjIxEaGgofv/9d7z88stIT0/H119/LWG0zpGWloa4uDiUlZVBrVZj8+bN6N69O44dO9bmzgl7YwG0rXNi/fr1OHr0KA4dOmTV58zfFUzCWpHRo0eb/9yrVy8MGDAAkZGR2LhxI6ZOnSphZNScTJgwwfzn2NhY9OrVC506dcLevXsxfPhwCSNznhkzZuD48eNtZo5kbeyNxdNPP23+c2xsLEJCQjB8+HBkZmaiU6dOTR2mU8XExODYsWMoLCzEl19+iaSkJKSkpEgdliTsjUX37t3bzDmRlZWF5557Djt27ICbm1uTfjZvR7ZiPj4+6NKlCzIyMqQORVLBwcEAYLWSJTc319zXlnXs2BHt2rVrtefJzJkzsXXrVuzZswfh4eHm9uDgYBiNRhQUFFhs35rPC3tjYcuAAQMAoFWeF0qlEtHR0ejTpw8WLlyI3r174x//+EebPCfsjYUtrfWcOHLkCPLy8nDbbbfB1dUVrq6uSElJwbJly+Dq6oqgoCCnnRdMwloxnU6HzMxMhISESB2KpKKiohAcHIxdu3aZ24qKinDw4EGLuQ9t1cWLF3H16tVWd54IITBz5kxs3rwZu3fvRlRUlEV/nz59oFAoLM6L9PR0XLhwodWdF3WNhS3Hjh0DgFZ3XthiMplgMBja1DlhT/VY2NJaz4nhw4cjLS0Nx44dM7/69u2LRx55xPxnZ50XvB3Zirzwwgu4//77ERkZiezsbMybNw8uLi6YOHGi1KE5nU6ns/jX2blz53Ds2DH4+fmhffv2mD17Nv72t7+hc+fOiIqKwmuvvYbQ0FCMGTNGuqCdpLax8PPzw4IFC5CYmIjg4GBkZmbipZdeQnR0NEaOHClh1I43Y8YMrFu3Dlu2bIGXl5d57oa3tzfc3d3h7e2NqVOnYu7cufDz84NGo8GsWbMQFxeH22+/XeLoHauuscjMzMS6detwzz33wN/fH7///jvmzJmDwYMHWy3Vb+mSk5MxevRotG/fHsXFxVi3bh327t2L7du3t6lzAqh9LNrSOeHl5WUxPxIAPD094e/vb2532nlxU2srqVkZP368CAkJEUqlUoSFhYnx48eLjIwMqcNqEnv27BEArF5JSUlCiKoyFa+99poICgoSKpVKDB8+XKSnp0sbtJPUNhYlJSVixIgRIiAgQCgUChEZGSmeeuopcfnyZanDdjhbYwBArF692rxNaWmpmD59uvD19RUeHh7iwQcfFDk5OdIF7SR1jcWFCxfE4MGDhZ+fn1CpVCI6Olq8+OKLorCwUNrAneCJJ54QkZGRQqlUioCAADF8+HDx008/mfvbyjkhRO1j0ZbOCVtqludw1nkhE0KIm0vjiIiIiKihOCeMiIiISAJMwoiIiIgkwCSMiIiISAJMwoiIiIgkwCSMiIiISAJMwoiIiIgkwCSMiIiISAJMwoiIiIgkwCSMiJqlyZMnN8ljpTp06IClS5c6/XOas7vuuguzZ8+WOgyiNodJGBHV2+TJkyGTySCTyaBUKhEdHY033ngDFRUVUodWpzVr1sDHx8eq/dChQ3j66aed/vkfffQRevfuDbVaDR8fH9x6661YuHCh0z+XiJovPsCbiBpk1KhRWL16NQwGA3744QfMmDEDCoUCycnJVtsajUYolUoJoqy/gIAAp3/GJ598gtmzZ2PZsmUYMmQIDAYDfv/9dxw/ftzpn01EzRevhBFRg6hUKgQHByMyMhLTpk1DfHw8vv32WwDXbyG+9dZbCA0NRUxMDAAgLS0Nw4YNg7u7O/z9/fH0009Dp9OZj1lZWYm5c+fCx8cH/v7+eOmll1Dzsba2bhvecsstmD9/vvl9QUEB/vKXvyAoKAhubm7o2bMntm7dir1792LKlCkoLCw0X8mr3q/mcS9cuICEhASo1WpoNBqMGzcOubm55v758+fjlltuwaeffooOHTrA29sbEyZMQHFxsd0x+/bbbzFu3DhMnToV0dHR6NGjByZOnIi33nrLvE312C1YsAABAQHQaDR45plnYDQazduYTCYsXLgQUVFRcHd3R+/evfHll19afNbx48cxevRoqNVqBAUF4bHHHsOVK1fM/Xq9Ho8//jjUajVCQkLw/vvv242biJyLSRgR3RR3d3eLRGHXrl1IT0/Hjh07sHXrVuj1eowcORK+vr44dOgQNm3ahJ07d2LmzJnmfd5//32sWbMGn3zyCfbt24dr165h8+bNDYrDZDJh9OjR2L9/Pz777DOcPHkSixYtgouLCwYOHIilS5dCo9EgJycHOTk5eOGFF2weIyEhAdeuXUNKSgp27NiBP/74A+PHj7fYLjMzE9988w22bt2KrVu3IiUlBYsWLbIbW3BwMFJTU3H+/Plav8OuXbtw6tQp7N27F1988QW+/vprLFiwwNy/cOFC/Oc//8GqVatw4sQJzJkzB48++ihSUlIAVCWhw4YNw6233orDhw9j27ZtyM3Nxbhx48zHePHFF5GSkoItW7bgp59+wt69e3H06NF6jTEROZggIqqnpKQkkZCQIIQQwmQyiR07dgiVSiVeeOEFc39QUJAwGAzmff71r38JX19fodPpzG3ff/+9kMvl4vLly0IIIUJCQsQ777xj7i8vLxfh4eHmzxJCiMjISLFkyRKLeHr37i3mzZsnhBBi+/btQi6Xi/T0dJuxr169Wnh7e1u133jcn376Sbi4uIgLFy6Y+0+cOCEAiF9//VUIIcS8efOEh4eHKCoqMm/z4osvigEDBtj8XCGEyM7OFrfffrsAILp06SKSkpLEhg0bRGVlpXmbpKQk4efnJ/R6vblt5cqVQq1Wi8rKSlFWViY8PDzEL7/8YnHsqVOniokTJwohhHjzzTfFiBEjLPqzsrIEAJGeni6Ki4uFUqkUGzduNPdfvXpVuLu7i+eee85u/ETkHJwTRkQNsnXrVqjVapSXl8NkMmHSpEkWtwRjY2Mt5oGdOnUKvXv3hqenp7ntjjvugMlkQnp6Otzc3JCTk4MBAwaY+11dXdG3b1+rW5K1OXbsGMLDw9GlS5dGf7dTp04hIiICERER5rbu3bvDx8cHp06dQr9+/QBU3cL08vIybxMSEoK8vDy7xw0JCcGBAwdw/Phx/Pzzz/jll1+QlJSEf//739i2bRvk8qqbEr1794aHh4d5v7i4OOh0OmRlZUGn06GkpAR33323xbGNRiNuvfVWAMBvv/2GPXv2QK1WW8WQmZmJ0tJSGI1Gi7H28/Mz3zYmoqbFJIyIGmTo0KFYuXIllEolQkND4epq+WvkxmTLkeRyuVVSVl5ebv6zu7u7Uz7XFoVCYfFeJpPBZDLVuV/Pnj3Rs2dPTJ8+Hc888wzuvPNOpKSkYOjQoXXuWz2H7vvvv0dYWJhFn0qlMm9z//334+2337baPyQkBBkZGXV+DhE1Hc4JI6IG8fT0RHR0NNq3b2+VgNnSrVs3/Pbbb9Dr9ea2/fv3Qy6XIyYmBt7e3ggJCcHBgwfN/RUVFThy5IjFcQICApCTk2N+X1RUhHPnzpnf9+rVCxcvXsSZM2dsxqFUKlFZWVlnrFlZWcjKyjK3nTx5EgUFBejevXud37Uhqo9347j89ttvKC0tNb9PTU2FWq1GREQEunfvDpVKhQsXLiA6OtriVX3l7rbbbsOJEyfQoUMHq208PT3RqVMnKBQKi7HWarV2x4yInItJGBE51SOPPAI3NzckJSXh+PHj2LNnD2bNmoXHHnsMQUFBAIDnnnsOixYtwjfffIPTp09j+vTpKCgosDjOsGHD8Omnn+K///0v0tLSkJSUBBcXF3P/kCFDMHjwYCQmJmLHjh04d+4cfvzxR2zbtg1A1S1EnU6HXbt24cqVKygpKbGKNT4+HrGxsXjkkUdw9OhR/Prrr3j88ccxZMgQ9O3bt9FjMG3aNLz55pvYv38/zp8/j9TUVDz++OMICAhAXFyceTuj0YipU6fi5MmT+OGHHzBv3jzMnDkTcrkcXl5eeOGFFzBnzhysXbsWmZmZOHr0KP75z39i7dq1AIAZM2bg2rVrmDhxIg4dOoTMzExs374dU6ZMQWVlJdRqNaZOnYoXX3wRu3fvxvHjxzF58mTz7VAialr8P4+InMrDwwPbt2/HtWvX0K9fPzz00EMYPnw4PvjgA/M2zz//PB577DEkJSUhLi4OXl5eePDBBy2Ok5ycjCFDhuC+++7DvffeizFjxqBTp04W23z11Vfo168fJk6ciO7du+Oll14yX/0aOHAgnnnmGYwfPx4BAQF45513rGKVyWTYsmULfH19MXjwYMTHx6Njx47YsGHDTY1BfHw8UlNT8fDDD6NLly5ITEyEm5sbdu3aBX9/f/N2w4cPR+fOnTF48GCMHz8eDzzwgMV8uzfffBOvvfYaFi5ciG7dumHUqFH4/vvvERUVBQAIDQ3F/v37UVlZiREjRiA2NhazZ8+Gj4+POdF69913ceedd+L+++9HfHw8Bg0ahD59+tzU9yOixpGJhsx8JSIip5g8eTIKCgrwzTffSB0KETURXgkjIiIikgCTMCIiIiIJ8HYkERERkQR4JYyIiIhIAkzCiIiIiCTAJIyIiIhIAkzCiIiIiCTAJIyIiIhIAkzCiIiIiCTAJIyIiIhIAkzCiIiIiCTw/wDXywBc0daNyQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 700x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "simulated_dataset[\"preds\"] = pps[\"obs\"].mean(axis=0)\n",
    "\n",
    "plt.figure(figsize=(7, 3))\n",
    "sns.histplot(data=simulated_dataset, x=\"Production_Speed\", binwidth=1, label=\"obs\")\n",
    "sns.histplot(data=simulated_dataset, x=\"preds\", binwidth=1, label=\"pps\")\n",
    "plt.xlabel(\"Production Speed\")\n",
    "plt.title(\"Posterior Predictive Distribution\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([13.455311, 10.143251], dtype=float32)"
      ]
     },
     "execution_count": 250,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "config_idx_new = jnp.array([0, 1], dtype=jnp.int32)\n",
    "\n",
    "# Correctly use a value of 0.0 for the third process parameter of product group 0\n",
    "X_test = jnp.array([\n",
    "    [3.0, 2.0, 0.0],\n",
    "    [3.0, 2.0, 2.0]\n",
    "])\n",
    "\n",
    "rng_key, rng_subkey = random.split(key=rng_key)\n",
    "pps_new = predictive(\n",
    "    rng_subkey,\n",
    "    parameter_mask,\n",
    "    config_idx_new,\n",
    "    X_test,\n",
    "    None\n",
    ")\n",
    "\n",
    "pps_new[\"obs\"].mean(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Average production speed for product group 0 is indeed higher than product group 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([20.98238 , 13.135274], dtype=float32)"
      ]
     },
     "execution_count": 253,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "config_idx_new = jnp.array([0, 1], dtype=jnp.int32)\n",
    "\n",
    "# Input a value of 2.0 for the third process parameter of product group 0\n",
    "# to see whether or not predictions are \"impacted\"\n",
    "X_test = jnp.array([\n",
    "    [6.0, 2.0, 2.0],\n",
    "    [6.0, 2.0, 2.0]\n",
    "])\n",
    "\n",
    "rng_key, rng_subkey = random.split(key=rng_key)\n",
    "pps_new = predictive(\n",
    "    rng_subkey,\n",
    "    parameter_mask,\n",
    "    config_idx_new,\n",
    "    X_test,\n",
    "    None\n",
    ")\n",
    "\n",
    "pps_new[\"obs\"].mean(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even with an incorrectly inputed value for the third process parameter, this value is effectively ignored resulting in an higher average production speed for product group 0 than product group 1."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".komax",
   "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.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}