COGMOD-HWI/HW04/HW4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2d616210",
   "metadata": {},
   "outputs": [],
   "source": [
    "# DATA\n",
    "\n",
    "import pandas as pd\n",
    "import stan\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import arviz as az\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "\n",
    "# stan problems\n",
    "import nest_asyncio\n",
    "nest_asyncio.apply()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c5cdedfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "url = \"https://raw.githubusercontent.com/stedy/Machine-Learning-with-R-datasets/master/insurance.csv\"\n",
    "data = pd.read_csv(url)\n",
    "\n",
    "# god is dead, force a normal\n",
    "data['smoker'] = data['smoker'].map({'yes': 1, 'no': 0})\n",
    "\n",
    "X = data[['bmi', 'age', 'children', 'smoker']].values\n",
    "y = data['charges'].values\n",
    "\n",
    "\n",
    "# Standardize predictors (mean=0, sd=1)\n",
    "X_mean = X.mean(axis=0)\n",
    "X_std = X.std(axis=0)\n",
    "y_mean = y.mean()\n",
    "y_std = y.std()\n",
    "y = (y - y_mean) / y_std\n",
    "X_std = X.std(axis=0)\n",
    "continuous_vars = data[['bmi', 'age', 'children']]\n",
    "continuous_vars_standardized = (continuous_vars - continuous_vars.mean()) / continuous_vars.std()\n",
    "\n",
    "X_standardized = np.hstack([continuous_vars_standardized, data[['smoker']].values])\n",
    "\n",
    "# Split into train/test\n",
    "X_train, X_test, y_train, y_test = train_test_split(X_standardized, y, test_size=0.2, random_state=42)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1a832a9e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Building...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Building: found in cache, done."
     ]
    }
   ],
   "source": [
    "stan_data = {\n",
    "  'N': X_train.shape[0],\n",
    "  'M': X_train.shape[1],\n",
    "  'X': X_train,\n",
    "  'y': y_train,\n",
    "  'N_test': X_test.shape[0],\n",
    "  'X_test': X_test\n",
    "}\n",
    "\n",
    "with open('multiple_regression.stan', 'r') as f:\n",
    "    stan_code = f.read()\n",
    "    \n",
    "posterior = stan.build(stan_code, data=stan_data, random_seed=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "def536a9",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling:   0%\n",
      "Sampling:  25% (3000/12000)\n",
      "Sampling:  50% (6000/12000)\n",
      "Sampling:  75% (9000/12000)\n",
      "Sampling: 100% (12000/12000)\n",
      "Sampling: 100% (12000/12000), done.\n",
      "Messages received during sampling:\n",
      "  Gradient evaluation took 0.00019 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 1.9 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:\n",
      "  Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/httpstan_gquumxii/model_rezrxovk.stan', line 21, column 2 to column 38)\n",
      "  If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,\n",
      "  but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.\n",
      "  Gradient evaluation took 0.000184 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 1.84 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Gradient evaluation took 0.000168 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 1.68 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Gradient evaluation took 0.000129 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 1.29 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:\n",
      "  Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/httpstan_gquumxii/model_rezrxovk.stan', line 21, column 2 to column 38)\n",
      "  If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,\n",
      "  but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n",
      "alpha   -0.396  0.017  -0.430   -0.364        0.0      0.0    9346.0   \n",
      "beta[0]  0.165  0.016   0.136    0.195        0.0      0.0   10163.0   \n",
      "beta[1]  0.298  0.016   0.269    0.328        0.0      0.0   10522.0   \n",
      "beta[2]  0.042  0.015   0.014    0.072        0.0      0.0   10397.0   \n",
      "beta[3]  1.951  0.039   1.875    2.020        0.0      0.0    9719.0   \n",
      "sigma    0.507  0.011   0.486    0.527        0.0      0.0    9713.0   \n",
      "\n",
      "         ess_tail  r_hat  \n",
      "alpha      6429.0    1.0  \n",
      "beta[0]    6037.0    1.0  \n",
      "beta[1]    5893.0    1.0  \n",
      "beta[2]    5565.0    1.0  \n",
      "beta[3]    6633.0    1.0  \n",
      "sigma      6303.0    1.0  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Posterior Distributions of Beta Coefficients')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGZCAYAAABxI8CQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPrhJREFUeJzt3XlclOX+//H3gAyrigbmUoK4ZKJlLmm550JmpYW5VKb2Fc0W86SlVh6wc4xj6dGyjpWmcigtl7RTqWgpaqappbZZoSm2uGC4Ysoy1+8Pf4yMLAIiAzev5+PBg5l77rnvz8w9y3uu+7qv22aMMQIAACjnPNxdAAAAQEkg1AAAAEsg1AAAAEsg1AAAAEsg1AAAAEsg1AAAAEsg1AAAAEsg1AAAAEsg1AAAAEsg1AAFmD9/vmw2m/bv3+/uUvI0ZMgQhYaGlsq6QkNDNWTIEOf17Odm+/btpbL+zp07q3PnzqWyritp27ZtuvXWW+Xv7y+bzaadO3e6uyRLyMzM1DPPPKNrr71WHh4e6tOnjyTp9OnTGjZsmGrWrCmbzabRo0dr//79stlsmj9/fpHWUdY/D0CoQQnLftNn//n4+KhRo0Z6/PHHdfjw4RJf35kzZxQTE6PExMQSX3Zpi4mJcXnu/Pz8VLduXd11112aN2+ezp07VyLr+eGHHxQTE1MmP5jLcm0lISMjQ/fdd59SU1M1ffp0xcfHKyQkJM95ExMTXV4PNptN1atXV9u2bfXuu+8Wu4YFCxZoxowZxb7/pSxbtkw9e/ZUUFCQ7Ha7ateurX79+mnt2rVXbJ2SNHfuXL388svq27ev4uLi9Le//U2S9OKLL2r+/PkaOXKk4uPjNWjQoCtax+Wy0meaWxigBM2bN89IMi+88IKJj483s2fPNoMHDzYeHh6mXr16Ji0trUTXl5KSYiSZ6OjoEl1utszMTPPXX38Zh8NxRZafU3R0tJFkZs2aZeLj482cOXPMpEmTzK233mokmRtuuMEcOHDA5T7p6enm7NmzRVrP4sWLjSSzbt26It3v7NmzJj093Xk9e1tv27atSMspbm3nzp0z586dK7F1ucPu3buNJDN79uxLzrtu3TojyYwaNcrEx8eb+Ph4M2PGDHPLLbcYSea1114rVg29evUyISEhxbpvQRwOhxkyZIiRZG666SYzefJk8/bbb5t//vOfpmXLlkaS2bRpU4mvN1v//v1NnTp1ck1v06aNadeuXa5a//rrL5OZmVmkdZTG58GV/kyzukpuSVKwvJ49e6pVq1aSpGHDhumqq67Sv//9b3344YcaOHCgm6u7tLS0NPn7+8vT01Oenp4lttwzZ87Iz8+vwHn69u2roKAg5/W///3vevfdd/XQQw/pvvvu05YtW5y3eXl5lVhteTHG6OzZs/L19ZW3t/cVXdel2O12t66/JBw5ckSSFBgYWOj7dOjQQX379nVeHzlypMLCwrRgwQI99thjJV1isU2bNk3z58/X6NGj9e9//1s2m81523PPPaf4+HhVqnTlvnKOHDmS5/N65MgRNWnSxGVadityUZX05wGuAHenKlhLfr/eP/74YyPJTJ482RhjTEZGhnnhhRdMWFiYsdvtJiQkxEyYMCFXq8O2bdtMjx49zFVXXWV8fHxMaGioGTp0qDHGmH379hlJuf5y/sLZvXu3iYyMNNWqVTPe3t6mZcuW5sMPP8yz5sTERDNy5EgTHBxsAgMDXW7bt2+fy31ef/1106RJE2O3202tWrXMo48+ao4dO+YyT6dOnUx4eLjZvn276dChg/H19TVPPvlkvs9ddktNSkpKnrcPHz7cSDKrV692Ths8eHCuX90LFy40LVq0MAEBAaZy5cqmadOmZsaMGS6P5+K/7JaRkJAQ06tXL7Nq1SrTsmVL4+3tbaZPn+68bfDgwbmet/Xr15vhw4eb6tWrm8qVK5tBgwaZ1NRUl5ou3i7Zci7zUrV16tTJdOrUyeX+hw8fNg8//LCpUaOG8fb2NjfccIOZP3++yzzZr5OXX37ZvPnmm87XXKtWrczWrVtd5j148KAZMmSIqVOnjrHb7aZmzZrm7rvvzrX98/LZZ5+Z9u3bGz8/P1O1alVz9913mx9++MF5++DBg3M9tosfT07ZLTWLFy/OdVvTpk1Nx44dc02Pj483LVq0MD4+PqZatWqmf//+Lq17nTp1ylVD9uvn3LlzZuLEiaZFixamSpUqxs/Pz7Rv396sXbv2ko/9zJkzpnr16qZx48aFbv3Yu3ev6du3r6lWrZrx9fU1bdq0MR9//HGu+c6ePWv+/ve/m/r16xu73W6uueYa8/TTTzs/K/L7HMh+/i7+27dvn/M+8+bNc1nX7t27zX333WeCgoKMj4+PadSokXn22Wedt+f3ebBixQrntg8ICDB33HGH+e6771zmGTx4sPH39ze//fab6d27t/H39zdBQUFmzJgxzufsUp9pl/P6rChoqUGp2Lt3ryTpqquuknS+9SYuLk59+/bVmDFj9OWXXyo2Nla7d+/WsmXLJJ3/hdWjRw8FBwdr/PjxCgwM1P79+/XBBx9IkoKDgzVr1iyNHDlS99xzj+69915J0g033CBJ+v7779WuXTvVqVNH48ePl7+/vxYtWqQ+ffpo6dKluueee1xqfPTRRxUcHKy///3vSktLy/exxMTEaNKkSerWrZtGjhypn376SbNmzdK2bdu0adMml9aTP//8Uz179tSAAQP04IMP6uqrry72czho0CC99dZbWr16tbp3757nPGvWrNHAgQPVtWtXTZkyRZK0e/dubdq0SU8++aQ6duyoUaNG6dVXX9Wzzz6r66+/XpKc/yXpp59+0sCBAzVixAhFRUXpuuuuK7Cuxx9/XIGBgYqJiXE+F8nJyc4+IYVVmNpy+uuvv9S5c2ft2bNHjz/+uOrVq6fFixdryJAhOn78uJ588kmX+RcsWKBTp05pxIgRstlseumll3Tvvffql19+cW6zyMhIff/993riiScUGhqqI0eOaM2aNTpw4ECBHbI//fRT9ezZU2FhYYqJidFff/2lmTNnql27dvr6668VGhqqESNGqE6dOnrxxRc1atQotW7dulCvh1OnTuno0aOSpNTUVC1YsEDfffed3n77bZf5Jk+erIkTJ6pfv34aNmyYUlJSNHPmTHXs2FE7duxQYGCgnnvuOZ04cUK//fabpk+fLkkKCAiQJJ08eVJz5szRwIEDFRUVpVOnTuntt99WRESEtm7dqubNm+db4+eff67U1FSNHj26UC0Zhw8f1q233qozZ85o1KhRuuqqqxQXF6e7775bS5Yscb43HQ6H7r77bn3++ecaPny4rr/+en377beaPn26fv75Zy1fvlzBwcGKj4/X5MmTdfr0acXGxko6/7qJj4/X3/72N11zzTUaM2aMpPOfGykpKblq+uabb9ShQwd5eXlp+PDhCg0N1d69e/XRRx9p8uTJ+T6W+Ph4DR48WBEREZoyZYrOnDmjWbNmqX379tqxY4fL6yYrK0sRERFq06aNpk6dqk8//VTTpk1T/fr1NXLkyEt+phX39VmhuDtVwVqyf8l8+umnJiUlxfz666/mvffeM1dddZXx9fU1v/32m9m5c6eRZIYNG+Zy37FjxxpJzl+Gy5Ytu2SfjYL2P3ft2tU0a9bMpfXH4XCYW2+91TRs2DBXze3bt8/1K/PiX2ZHjhwxdrvd9OjRw2RlZTnne+2114wkM3fuXOe07F/Fb7zxxqWfOHPplppjx44ZSeaee+5xTru4pebJJ580VapUKfDXckH9VkJCQowks2rVqjxvy6ulpmXLli59bV566SUjyaVFLL9tdPEyC6rt4paaGTNmGEnmnXfecU5LT083t9xyiwkICDAnT540xlz49XvVVVe5tCB9+OGHRpL56KOPjDEXnt+XX34517ovpXnz5qZGjRrmzz//dE7btWuX8fDwMA899JBzWkGtLxfLr6XBw8PD2eKZbf/+/cbT0zPX9G+//dZUqlTJZXp+fWoyMzNz9Vk6duyYufrqq83DDz9cYK2vvPKKkWSWLVt2ycdljDGjR482kszGjRud006dOmXq1atnQkNDne+t+Ph44+Hh4TKfMca88cYbufroZLeMXiy79TGnvFpqOnbsaCpXrmySk5Nd5s3Zf+biz4NTp06ZwMBAExUV5XKfQ4cOmapVq7pMz26pe+GFF1zmvemmm0zLli2d1/P7TLuc12dFwtFPuCK6deum4OBgXXvttRowYIACAgK0bNky1alTRytWrJAkPfXUUy73yf4l9cknn0i60O/g448/VkZGRpHWn5qaqrVr16pfv37OX7pHjx7Vn3/+qYiICCUlJen33393uU9UVNQlf2V++umnSk9P1+jRo+XhceHtExUVpSpVqjhrz+bt7a2hQ4cWqfb8ZP+iPnXqVL7zBAYGKi0tTWvWrCn2eurVq6eIiIhCzz98+HCX1qmRI0eqUqVKzu18paxYsUI1a9Z06aPl5eWlUaNG6fTp01q/fr3L/P3791e1atWc1zt06CBJ+uWXXyRJvr6+stvtSkxM1LFjxwpdx8GDB7Vz504NGTJE1atXd06/4YYb1L1798t+Hv7+979rzZo1WrNmjd5//30NHDhQzz33nF555RXnPB988IEcDof69evnfK0fPXpUNWvWVMOGDbVu3bpLrsfT09PZb8nhcCg1NVWZmZlq1aqVvv766wLve/LkSUlS5cqVC/WYVqxYoZtvvlnt27d3TgsICNDw4cO1f/9+/fDDD5KkxYsX6/rrr1fjxo1dHtdtt90mSYV6XIWRkpKiDRs26OGHH1bdunVdbiuotXHNmjU6fvy4Bg4c6FKfp6en2rRpk2d9jzzyiMv1Dh06OF+DBSnu67OiYfcTrojXX39djRo1UqVKlXT11Vfruuuuc4aA5ORkeXh4qEGDBi73qVmzpgIDA5WcnCxJ6tSpkyIjIzVp0iRNnz5dnTt3Vp8+fXT//fdfstPqnj17ZIzRxIkTNXHixDznOXLkiOrUqeO8Xq9evUs+ruzaLt4lY7fbFRYW5rw9W506dUqsg+vp06clFfzF8eijj2rRokXq2bOn6tSpox49eqhfv366/fbbC72ewjwPOTVs2NDlekBAgGrVqnXFD8tOTk5Ww4YNXcKldGF31cXb4uIvq+yAk/0F4e3trSlTpmjMmDG6+uqr1bZtW91555166KGHVLNmzQLrkHK/JrJrSUhIcHY8L45mzZqpW7duzuv9+vXTiRMnNH78eN1///0KDg5WUlKSjDG5tkW2wnYoj4uL07Rp0/Tjjz+6/JC41GuiSpUqkgoO3DklJyerTZs2uabn3HZNmzZVUlKSdu/ereDg4DyXk93x+nJlh4qmTZsW6X5JSUmS5AxZF8t+XrL5+PjkeizVqlUrVEgp7uuzoiHU4Iq4+eabnUc/5edS/S1sNpuWLFmiLVu26KOPPlJCQoIefvhhTZs2TVu2bHG2XOTF4XBIksaOHZtvq8PFocrX17fAeoqjJJf53XffScpdd041atTQzp07lZCQoJUrV2rlypWaN2+eHnroIcXFxRVqPVfiechPVlZWqa0rv1Y4Y4zz8ujRo3XXXXdp+fLlSkhI0MSJExUbG6u1a9fqpptuKq1SL6lr1676+OOPtXXrVvXq1UsOh0M2m00rV67M83EW9F7J9s4772jIkCHq06ePnn76adWoUUOenp6KjY119onLT+PGjSVJ3377rXPQu5LgcDjUrFkz/fvf/87z9muvvbbE1lUc2Z8z8fHxeQaLi4/2utwjp8rL69OdCDUodSEhIXI4HEpKSnLpBHr48GEdP34812Bkbdu2Vdu2bTV58mQtWLBADzzwgN577z0NGzYs32AUFhYm6fwv1Jy/ckuidul8Z9rsdUhSenq69u3bV6Lrulh8fLwkXXLXkN1u11133aW77rpLDodDjz76qN58801NnDhRDRo0KFLn3cJISkpSly5dnNdPnz6tgwcP6o477nBOq1atmo4fP+5yv/T0dB08eNBlWlFqCwkJ0TfffCOHw+HSWvPjjz86by+O+vXra8yYMRozZoySkpLUvHlzTZs2Te+8806+dUjnXxMX+/HHHxUUFFTsVpr8ZGZmSrrQele/fn0ZY1SvXj01atSowPvm9xwvWbJEYWFh+uCDD1zmiY6OvmQ97du3V7Vq1bRw4UI9++yzl/zyDgkJyff5yr5dOv+4du3apa5du5b46zan7Pdy9g+Hwqpfv76k8z8mSuq9f6nHWdTXZ0VDnxqUuuwvu4tHNc3+NdarVy9J53cL5PwVLcl5BEb26LrZY75c/IVZo0YNde7cWW+++WauL05JeR79UBjdunWT3W7Xq6++6lLb22+/rRMnTjhrL2kLFizQnDlzdMstt6hr1675zvfnn3+6XPfw8HAeOZH9nGV/wV78nBXXW2+95bKrYtasWcrMzFTPnj2d0+rXr68NGzbkut/FLTVFqe2OO+7QoUOH9P777zunZWZmaubMmQoICFCnTp2K9DjOnDmjs2fPukyrX7++KleuXOBozrVq1VLz5s0VFxfnUvd3332n1atXu4S7kvLxxx9Lkm688UZJ0r333itPT09NmjQp13vGGOPyuvD399eJEydyLTM7iOS8/5dffqnNmzdfsh4/Pz+NGzdOu3fv1rhx43LVIJ1vCdq6dauk89tu69atLstOS0vTW2+9pdDQUOe4Mv369dPvv/+u2bNn51reX3/9VeBRikURHBysjh07au7cuTpw4IDLbXk9lmwRERGqUqWKXnzxxTz7/RXncya/z7Tivj4rGlpqUOpuvPFGDR48WG+99ZaOHz+uTp06aevWrYqLi1OfPn2cv/rj4uL0n//8R/fcc4/q16+vU6dOafbs2apSpYrzi8LX11dNmjTR+++/r0aNGql69epq2rSpmjZtqtdff13t27dXs2bNFBUVpbCwMB0+fFibN2/Wb7/9pl27dhW59uDgYE2YMEGTJk3S7bffrrvvvls//fST/vOf/6h169Z68MEHL/v5WbJkiQICApSenq7ff/9dCQkJ2rRpk2688UYtXry4wPsOGzZMqampuu2223TNNdcoOTlZM2fOVPPmzZ2tYs2bN5enp6emTJmiEydOyNvbW7fddptq1KhRrHrT09PVtWtX9evXz/lctG/fXnfffbdLXY888ogiIyPVvXt37dq1SwkJCS6DDBa1tuHDh+vNN9/UkCFD9NVXXyk0NFRLlizRpk2bNGPGjEJ3Ws32888/Ox9HkyZNVKlSJS1btkyHDx/WgAEDCrzvyy+/rJ49e+qWW27R//3f/zkP6a5atapiYmKKVMfFNm7c6PwyS01N1f/+9z+tX79eAwYMcO72qV+/vv75z39qwoQJ2r9/v/r06aPKlStr3759WrZsmYYPH66xY8dKklq2bKn3339fTz31lFq3bq2AgADddddduvPOO/XBBx/onnvuUa9evbRv3z698cYbatKkibNFqCBPP/20vv/+e02bNk3r1q1T3759VbNmTR06dEjLly/X1q1b9cUXX0iSxo8fr4ULF6pnz54aNWqUqlevrri4OO3bt09Lly51trwNGjRIixYt0iOPPKJ169apXbt2ysrK0o8//qhFixYpISHhkru5C+vVV19V+/bt1aJFCw0fPlz16tXT/v379cknn+R7fq4qVapo1qxZGjRokFq0aKEBAwYoODhYBw4c0CeffKJ27drptddeK1Id+X2mZWZmFvv1WaG46agrWFRhh87PyMgwkyZNMvXq1TNeXl7m2muvzTX43tdff20GDhxo6tata7y9vU2NGjXMnXfeabZv3+6yrC+++MK0bNnS2O32XIdC7t271zz00EOmZs2axsvLy9SpU8fceeedZsmSJYWqOb/Btl577TXTuHFj4+XlZa6++mozcuTIfAffK6zsQ7qz/3x8fMw111xj7rzzTjN37tw8T4dw8SHdS5YsMT169DA1atQwdrvd1K1b14wYMcIcPHjQ5X6zZ882YWFhxtPTM8/B9/JyqcH3qlWrZgICAswDDzzgcmizMcZkZWWZcePGmaCgIOPn52ciIiLMnj17ci2zoNryG3xv6NChJigoyNjtdtOsWbNcA6rlHHzvYjlfL0ePHjWPPfaYady4sfH39zdVq1Y1bdq0MYsWLcrz+bjYp59+atq1a2d8fX1NlSpVzF133eUy+J4xl39It91uN40bNzaTJ092OYw+29KlS0379u2Nv7+/8ff3N40bNzaPPfaY+emnn5zznD592tx///0mMDDQZfA9h8NhXnzxRRMSEmK8vb3NTTfdZD7++OM8B3gsSPZrsHr16qZSpUqmVq1apn///iYxMdFlvuzB9wIDA42Pj4+5+eab8xx8Lz093UyZMsWEh4cbb29vU61aNdOyZUszadIkc+LECed8l3tItzHGfPfdd+aee+5x1nTdddeZiRMnOm/P7/Ng3bp1JiIiwlStWtX4+PiY+vXrmyFDhrh8VmUPvnex7Pd9Tnl9pl3u67OisBlTQNsaAABAOUGfGgAAYAmEGgAAYAmEGgAAYAmEGgAAYAmEGgAAYAmMU1NKjDGFPi8KAADIrXLlygWOukyoKSWnTp1S1apV3V0GAADl1okTJ3KdKDQnxqkpJSXRUpOenq5p06ZJksaMGVNiZ38GAKA8oKWmjLDZbAWmy8JIT0+Xj4+PpPPDcxNqAAC4gI7CAADAEgg1AADAEgg1AADAEgg1AADAEugoXI7YbDaFhIQ4LwMAgAs4pBsAAFgCu58AAIAlEGoAAIAl0KemHElPT9crr7wiSXryyScZfA8AgBwINeXMmTNn3F0CAABlErufAACAJRBqAACAJRBqAACAJRBqAACAJRBqAACAJXD0Uzlis9lUu3Zt52UAAHABp0kAAACWwO4nAABgCYQaAABgCfSpKUcyMjL0+uuvS5Iee+wxeXl5ubkiAADKDkJNOWKM0YkTJ5yXAQDABex+AgAAlkCoAQAAlkCoAQAAlkCoAQAAlkCoAQAAlsDRT+WIzWZTcHCw8zIAALiA0yQAAABLYPcTAACwBEINAACwBPrUlCMZGRmaPXu2JCkqKorTJAAAkAOhphwxxiglJcV5GQAAXMDuJwAAYAmEGgAAYAmEGgAAYAmEGgAAYAmEGgAAYAkc/VSO2Gw2Va1a1XkZAABcwGkSAACAJbD7CQAAWAKhBgAAWAJ9asqRjIwMzZ8/X5I0ZMgQTpMAAEAOhJpyxBijP/74w3kZAABcwO4nAABgCYQaAABgCYQaAABgCYQaAABgCYQaAABgCRz9VM74+fm5uwQAAMokTpMAAAAsgd1PAADAEgg1AADAEuhTU45kZGTo3XfflSQ98MADnCYBAIAcCDXliDFGycnJzssAAOACdj8BAABLINQAAABLINQAAABLINQAAABLINQAAABL4OincobDuAEAyBunSQAAAJbA7icAAGAJhBoAAGAJ9KkpRzIzM7Vo0SJJUr9+/VSpEpsPAIBsfCuWIw6HQ0lJSc7LAADgAnY/AQAASyDUAAAASyDUAAAASyDUAAAASyDUAAAASyDUAAAASyhyqElMTJTNZlNMTMwVKKdkxcTEyGazOf/Gjx9/WcsbP368y/JK+zmw2+2Kjo5WdHS07HZ7qa4bAICyrsyNUxMaGipJ2r9/f4ktc/DgwQoNDVX79u1z3Xby5EnFxMRo6dKlOnTokGrVqqX77rtP0dHRCggIcJm3W7du8vHx0f79+xUXF1di9QEAgMtX5kLNlTBkyBB17tw51/S0tDR16tRJO3fuVI8ePTRw4EDt2LFDU6dO1fr167Vhwwb5+Pg45+/WrZu6deumxMREQo2klB9SlLonVdUbVFdwk2B3lwMAqOAqRKjJz0svvaSdO3dq3Lhx+te//uWcPn78eE2ZMkXTp0/XhAkT3Fihq8zMTC1btkySdM8997jtNAlnjp7RoshFSt6Q7JwW1i1MkQsj5Rfk55aaAAC4rI7Cn3/+uTp37qzKlSsrMDBQkZGR2rNnT675jhw5or/97W9q0KCBvL29FRQUpMjISH333XfOefbv3y+bzabk5GQlJyfn2XclPT1dM2fOVEREhK699lp5e3urRo0auvfee7Vjx44i1W6M0Zw5cxQQEKCJEye63DZx4kQFBARozpw5RX9SriCHw6EffvhBP/zwQ6mfJiE9Ld35t7jfYh3+9rD6LuqrMQfHqO+ivjq065AW91/snAcAgNJW7J/6W7ZsUWxsrG6//XY98cQT+v7777Vs2TJt3LhRW7ZsUVhYmCRp79696ty5s3777Tf16NFDffr00ZEjR7R06VIlJCTos88+U5s2bRQYGKjo6GjNmDFDkjR69GjnurJ3HaWmpmr06NHq0KGD7rjjDlWrVk2//PKL/ve//2nlypXasGGDWrduXaj6k5KS9McffygiIkL+/v4ut/n7+6tdu3ZKSEjQr7/+qmuvvba4T5NlxAbEulzvu6ivwu8Ll6Tz/420pP8S53zRJrrUawQAVGzFDjUJCQl64403NGLECOe0N998U4888oiefPJJffTRR5Kkhx56SAcPHtSqVasUERHhnPf5559Xq1atFBUVpW+++UaBgYGKiYnR/PnzJSnPI4uqVaumAwcOqE6dOi7Tv//+e7Vt21bPPvus1qxZU6j6s08M2bBhwzxvb9iwoRISEpSUlESoyUNIhxDX6x1D8pkTAIDSUezdT40aNVJUVJTLtKioKDVs2FCffPKJUlJStGPHDn3xxRcaPHiwS6DJef9vv/3WZTdUQby9vXMFGkkKDw9Xly5dtGHDBmVkZBRqWSdOnJAkVa1aNc/bq1Sp4jJfRTfh9ARNOD1BUdvOb/Pkjckut2f3r4naFqUJp8tOPyQAQMVR7Jaadu3aycPDNRN5eHioXbt2SkpK0q5du5ytIYcPH86z5eXHH390/m/atGmh1rtz50699NJL+vzzz3Xo0KFcIebo0aOqVatWMR4RCmL3Pz8uTu1WtRXWLUwrHlshmfMtNMkbkrXi8RUK6x6m2q1qu7lSAEBFVexQc/XVVxc4/cSJE0pNTZUkffLJJ/rkk0/yXVZaWlqh1vnFF1/otttukyT16NFDDRs2VEBAgGw2m5YvX65du3bp3LlzhVpWdgtNfi0xJ0+edJkPF0QujNTS+5dqSf8lzmlh3cMUuSDSjVUBACq6Yoeaw4cPFzi9atWqzl04M2fO1OOPP17cVTlNnjxZ586d08aNG3MNpLdlyxbt2rWr0MvK7kuT3Zp0sUv1uanI/IL8NGj1IMapAQCUKcUONZs2bZLD4XDZBeVwOPTFF1/IZrPpxhtvdIaazZs3FzrUeHp6Kj0970OC9+7dq+rVq+cKNGfOnNHXX39dpPobNmyo2rVra9OmTUpLS3M5AiotLU2bNm1SvXr1ylQnYS8vL+e4OV5eXm6uRgpuEkyYAQCUGcXuKPzzzz9r9uzZLtNmz56tn3/+Wb169VJwcLBuvvlmtWnTRgsXLtT777+faxkOh0Pr1693mVa9enUdPXpUZ8+ezTV/SEiIjh07pu+//945LSsrS2PHjlVKSkqR6rfZbBo2bJhOnz6tf/zjHy63/eMf/9Dp06dzdYR2N5vNJrvdLrvdLpvN5u5yAAAoU4rdUhMREaFRo0ZpxYoVCg8P1/fff6+PPvpIQUFBeuWVV5zzLVy4UF26dNGAAQM0Y8YMtWjRQr6+vjpw4IA2b96slJQUlwBz2223afv27erZs6c6dOggu92ujh07qmPHjnriiSe0evVqtW/fXv369ZOPj48SExP1+++/q3PnzkpMTCzSY3jmmWf04YcfasqUKdqxY4datGihr7/+WqtXr1br1q1dxsoBAABlW7Fbatq2bavPPvtMJ06c0KuvvqrExET16dNHmzdvdg68J0n16tXTjh079Pzzz+v06dOaN2+e3nzzTe3cuVMdO3bUwoULXZY7ceJERUVF6aefftKLL76oiRMnau3atZKkO++8U0uWLFFYWJjeeecdLViwQI0bN9bWrVsVElL0cVL8/f21fv16jR49Wrt379a0adP0448/asyYMfrss8/k6+tb3KfnisjMzNTy5cu1fPlyZWZmurscAADKFJsxxri7iCslJiZGkyZN0rp16/I8oWVxJSYmqkuXLoqOjs7zUPUrJT09XbGx50fsnTBhgux2e6mtGwCAsu6yzv1UXnTp0kU2m03jx4+/rOWMHz9eNptNXbp0KaHKAABASbH0Wbovbp25+KipourWrZt8fHzyXT4AAHAfS+9+shp2PwEAkL8KsfsJAABYH6EGAABYAqEGAABYAn1qyhFjjM6cOSNJ8vPzY1RhAAByINQAAABLYPcTAACwBEuPU2M1mZmZSkhIkHT+3FuVKrH5AADIRktNOeJwOLR9+3Zt375dDofD3eUAAFCmEGoAAIAlEGoAAIAlEGoAAIAlEGoAAIAlEGoAAIAlEGoAAIAlMKJwOWKM0YkTJyRJVatW5TQJAADkQKgBAACWwO4nAABgCYyzX45kZWXps88+kyR17dpVnp6ebq4IAICyg5aaciQrK0ubN2/W5s2blZWV5e5yAAAoUwg1AADAEgg1AADAEgg1AADAEgg1AADAEgg1AADAEgg1AADAEhhRuBwxxiglJUWSFBwczGkSAADIgVADAAAsgd1PAADAEjhNQjmSlZWljRs3SpI6dOjAaRIAAMiBUFOOZGVlaf369ZKkW2+9lVADAEAO7H4CAACWQKgBAACWQKgBAACWQKgBAACWQKgBAACWQKgBAACWwIjC5YjD4dDBgwclSbVq1ZKHB5kUAIBshBoAAGAJ/NQHAACWwIjC5UhWVpa2bNkiSWrbti0jCgMAkAOhphzJysrSp59+Kklq3bo1oQYAgBzY/QQAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACyBEYXLEYfDoQMHDkiS6taty2kSAADIgVADAAAsgZ/6AADAEhhRuBzJysrSV199JUlq2bIlIwoDAJADoaYcycrK0sqVKyVJzZs3J9QAAJADu58AAIAlEGoAAIAlWDrUzJ8/Xzabzfk3YMCAYi9r1apVLsvq3LlzyRUKAAAuW4XoU9O7d281b95cTZs2dU7bu3ev4uPj9fXXX+urr77SH3/8oZCQEO3fvz/PZTRo0EDR0dGSpEmTJpVG2QAAoAgqRKjp06ePhgwZ4jJt48aNmjRpkjw9PXX99dfr0KFDBS6jQYMGiomJkUSocZeUH1KUuidV1RtUV3CTYHeXAwAoYypEqMlLx44dtXnzZt14443y9fWVj4+Pu0tCPs4cPaNFkYuUvCHZOS2sW5giF0bKL8jPjZUBAMqSChtqwsLCFBYW5u4yiqRSpUoaOHCg87KVpaelOy8v7rdYh789rL6L+iqkQ4iSNyZrxWMrtLj/Yg383/nnw+5vd1epAIAywtrfjBbj4eGhRo0aubuMUhEbEOtyve+ivgq/L1ySzv830pL+S5zzRZvoUq8RAFC2WProJ1hHSIcQ1+sdQ/KZEwBQUdFSU45kZWXp22+/lSQ1a9bM0iMKTzg9QZJ0dPdRzW49W8kbk50tNZKc/WuitkUp6Pogt9QIAChbCDXlSFZWlj788ENJUpMmTSwdarL7yNRuVVth3cK04rEVkjnfQpO8IVkrHl+hsO5hqt2qtpsrBQCUFYQalHmRCyO19P6lWtJ/iXNaWPcwRS6IdGNVAICyhlCDMs8vyE+DVg9inBoAQIEINSg3gpsEE2YAAPni6CcAAGAJFbal5ujRoxo7dqzzekZGho4ePepyOoWpU6cqKIgjawAAKA8qbKg5ffq04uLiXKalpaW5TIuJiSHUAABQTlTYUBMaGipjjLvLKJJKlSqpb9++zssAAOCCCtGnZujQobLZbBowYECxl7Fq1SrZbDbZbLYSrKxoPDw8FB4ervDwcHl4VIhNBwBAoVn6537z5s0VHX3hnEBNmzYt9rIaNGjgsqzQ0NDLKQ0AAJQwmylv+2AqMIfDod27d0uSrr/+elprAADIgW/FciQzM1NLlizRkiVLlJmZ6e5yAAAoUwg1AADAEgg1AADAEgg1AADAEgg1AADAEgg1AADAEgg1AADAEiw9+J7VeHp6qnfv3s7LAADgAgbfAwAAlsDuJwAAYAnsfipHHA6H9uzZI+n8uag4TQIAABfwrViOZGZmauHChVq4cCGnSQAA4CKEGgAAYAmEGgAAYAmEGgAAYAmEGgAAYAmEGgAAYAmEGgAAYAmMU1OOeHp6qmfPns7LAADgAk6TAAAALIHdTwAAwBLY/VSOOBwOHThwQJJUt25dTpMAAEAOfCuWI5mZmYqLi1NcXBynSQAA4CKEGgAAYAmEGgAAYAmEGgAAYAmEGgAAYAmEGgAAYAmEGgAAYAmMU1OOeHp6qlu3bs7LAADgAk6TAAAALIHdTwAAwBLY/VSOOBwOHTx4UJJUq1YtTpMAAEAOfCuWI5mZmZozZ47mzJnDaRIAALgIoQYAAFgCoQYAAFgCoQYAAFgCoQYAAFgCoQYAAFgCoQYAAFgC49SUI56enurUqZPzMgAAuIDTJAAAAEtg9xMAALAEdj+VI8YYpaSkSJKCg4Nls9ncXBEAAGUHLTXlSEZGhmbNmqVZs2YpIyPD3eUAAFCmEGoAAIAlEGoAAIAlEGoAAIAlEGoAAIAlEGoAAIAlEGoAAIAlME5NOeLp6albbrnFeRkAAFzAaRIAAIAlsPsJAABYgqVDzfz582Wz2Zx/AwYMKPayVq1a5bKszp07l1yhhWSM0fHjx3X8+HHRwAYAgKsK0aemd+/eat68uZo2bSrpfDhYtWqV/ve//2nTpk1KTk5WRkaGGjZsqP79++upp56Sj4+PyzIaNGig6OhoSdKkSZNK/TFI50+T8Morr0iSJkyYILvd7pY6AAAoiypEqOnTp4+GDBnivH7u3Dndcccd8vb2VufOnRUREaGzZ88qISFBzz33nJYvX67ExET5+fk579OgQQPFxMRIcl+oKayUH1KUuidV1RtUV3CTYHeXAwBAqagQoeZinp6e+uc//6lHH31U1apVc07PyMhQZGSkPvroI73++ut6+umn3Vhl0Z05ekaLIhcpeUOyc1pYtzBFLoyUX5BfAfcEAKD8s3Sfmvx4eXnpueeecwk02dMnTJggSVq/fr07Siuy9LR059/ifot1+NvD6ruor8YcHKO+i/rq0K5DWtx/sXMeAACsqkK21BTEy8tLklSpUvl4amIDYl2u913UV+H3hUvS+f9GWtJ/iXO+aBNd6jUCAFAaKmRLTUHmzp0rSerRo4ebKymekA4hrtc7huQzJwAA1lI+miNKycqVK/Xmm2/q+uuv1//93/+5u5xCmXD6/O6yo7uPanbr2UremOxsqZHk7F8TtS1KQdcHuaVGAABKA6Hm/9u2bZv69++vqlWravHixfL29nZ3Sbl4eHioVatWzsuSZPc/f1h37Va1FdYtTCseWyGZ8y00yRuSteLxFQrrHqbarWq7rW4AAEoDoUbS9u3b1aNHD3l4eCghIUHh4eGXvpMbVKpUSb169cr39siFkVp6/1It6b/EOS2se5giF0SWRnkAALhVhQ8127dvV/fu3eVwOLR69Wq1bt3a3SUVm1+QnwatHsQ4NQCACqlCh5rsQJOVlaWEhAS1adPG3SUVyBijM2fOSJL8/Pxks9nynC+4STBhBgBQ4VTYo5+++uorde/eXZmZmVq5cqVuueUWd5d0SRkZGZo6daqmTp2qjIwMd5cDAECZUiFbalJTU9W9e3cdP35ct99+u9asWaM1a9a4zBMYGKjRo0e7p0AAAFBkFTLUnDx5UseOHZN0/uzbq1atyjVPSEgIoQYAgHKkQoaa0NBQGWPcXQYAAChBFaJPzdChQ2Wz2TRgwIBiL2PVqlWy2Wz5ds4FAADuZemWmubNmys6+sK5jpo2bVrsZTVo0MBlWaGhoZdTGgAAKGE2w36YciM9PV2xsedPTDlhwgTZ7XY3VwQAQNlh6ZYaq/Hw8NCNN97ovAwAAC6gpQYAAFgCP/cBAIAlsPupHDHGOEcS9vLy4kgsAAByoKWmHMnIyFBsbKxiY2M5TQIAABch1AAAAEsg1AAAAEsg1AAAAEsg1AAAAEsg1AAAAEsg1AAAAEtgnJpyxMPDQ02aNHFeBgAAF3CaBAAAYAn83AcAAJZAqAEAAJZAn5pyJD09XbGxsZKkCRMmyG63u7kiAADKDlpqAACAJRBqAACAJRBqAACAJRBqAACAJRBqAACAJRBqAACAJXBIdzni4eGhhg0bOi8DAIALOE0CAACwBH7uAwAASyDUAAAAS6BPTTmSnp6uqVOnSpLGjh3LaRIAAMiBUFPOZGRkuLsEAADKJHY/AQAASyDUAAAASyDUAAAASyDUAAAASyDUAAAAS+Dop3LEZrMpJCTEeRkAAFzAaRIAAIAlsPsJAABYAqEGAABYAn1qypH09HS98sorkqQnn3yS0yQAAJADoaacOXPmjLtLAACgTGL3EwAAsARCDQAAsARCDQAAsARCDQAAsARCDQAAsASOfipHbDabateu7bwMAAAu4DQJAADAEtj9BAAALIFQAwAALIE+NeVIRkaGXn/9dUnSY489Ji8vLzdXBABA2UGoKUeMMTpx4oTzMgAAuIDdTwAAwBIINQAAwBIINQAAwBIINQAAwBIsHWrmz58vm83m/BswYECxl7Vq1SqXZXXu3LnkCgUAAJetQhz91Lt3bzVv3lxNmzZ1Tlu5cqXi4uK0c+dOHTp0SOnp6apbt67atWuncePGqVGjRi7LaNCggaKjoyVJkyZNKtX6s9lsNgUHBzsvAwCACyx9moT58+dr6NChmjdvnoYMGeJy2xNPPKGPPvpIbdq0Ue3ateXl5aXdu3dr5cqVqlSpklasWKHbbrstz+XabDZ16tRJiYmJV/5BAABQylJ+SFHqnlRVb1BdwU2C3V1OoVWIlpq8vPzyy5o5c2au6Z999pm6deumcePGadu2bW6oDACA0pOWkua8fObPM/pkxCdK3pDsnBbSMUS93uwlv6v8JEn+wf6lXmNhVdhQ4+Pjk+f0rl27qlq1atqzZ08pVwQAQMlJT0sv1HxTa0x1XrZ52ORd1Vt9F/VVSIcQJW9M1scjPtas8FkyjvM7diacnnDJZdr97cUr+jJV2FCTn82bN+vYsWNq3769u0vJJSMjQ7Nnz5YkRUVFcZoEAEC+YgNii3wf4zC68807FX5fuCSd/2+kJf2XFGm50Sa6yOsuCRU+1KxevVpffPGFzp07p6SkJH388ccKCgrS9OnT3V1aLsYYpaSkOC8DAFDSQjqEuF7vGJLPnGUPoWb1ak2bNs15vUGDBnrvvffUsmVLN1YFAMDlKcxuIil3y0vyxmRnS40kl/41RVmuO1T4UDN16lRNnTpVp0+f1g8//KAXXnhB7dq109y5c3X//fe7uzwAAIqlsP1axh4Z67y8uO9irXhshWTOt9Akb0jWisdXKKRTiO5bfF+RlusOFT7UZAsICNDNN9+s5cuXq1WrVho+fLi6d+/uHBcGAAArynk0U7+l/bT0/qUufWjCuocpckGk/IL83FFekRBqLlKpUiV16dJFu3bt0vbt29WzZ093lwQAQKnwC/LToNWDGKfGSv744w9J4ugiAECFFNwkuFyFmWyWPvdTQbZv357n9ISEBC1btkyBgYG65ZZbSrmqgtlsNlWtWlVVq1blNAkAAFykwrbUtG7dWk2bNtUNN9yga665Rmlpafrmm2+0ceNGeXl5ae7cufL3L1ujJnp5eWn06NHuLgMAgDKpwoaaF198UevWrdP69euVkpIiDw8P1a1bV8OHD9fo0aN1/fXXu7tEAABQBBU21EyYMEETJpTdY+0BAEDRVIg+NUOHDpXNZtOAAQOKvYxVq1bJZrO5tS9L9mkSZs+erYyMDLfVAQBAWWTplprmzZsrOvrC+SeaNm1a7GU1aNDAZVmhoaGXU1qxGGOcR2ZxmgQAAFzZDN+O5UZ6erpiY88PZz1hwgTZ7WV3VEcAAEpbhdj9BAAArI9QAwAALIFQAwAALIFQAwAALMHSRz9ZkZ9f2T9LKgAA7sDRTwAAwBLY/QQAACyBUAMAACyBPjXlSEZGht59911J0gMPPCAvLy83VwQAQNlBqClHjDFKTk52XgYAABew+wkAAFgCoQYAAFgCoQYAAFgCoQYAAFgCoQYAAFgCRz+VMxzGDQBA3jhNAgAAsAR2PwEAAEsg1AAAAEugT005kpmZqUWLFkmS+vXrp0qV2HwAAGTjW7EccTgcSkpKcl4GAAAXsPsJAABYAqEGAABYAqEGAABYAqEGAABYAqEGAABYAkc/lRJjjE6dOnVZy0hPT9fZs2clSSdPnpTdbi+J0gAAKBcqV64sm82W7+2cJqGUnDx5UlWrVnV3GQAAlFsnTpxQlSpV8r2dUFNKSqKlRjofjq699lr9+uuvBW5YuA/bqHxgO5V9bKOyr7S30aVaatj9VEpsNluJbvAqVarwJi/j2EblA9up7GMblX1lZRvRURgAAFgCoQYAAFgCoaac8fb2VnR0tLy9vd1dCvLBNiof2E5lH9uo7Ctr24iOwgAAwBJoqQEAAJZAqAEAAJZAqAEAAJZAqAEAAJZAqAEAAJZAqCkDtm3bpjvuuEOBgYHy9/dX27ZttWjRoiIt49y5c3rhhRfUsGFD+fj4qHbt2ho+fLiOHDlyhaquWC53G82fP182my3fv8TExCtXfAXwzjvvaMSIEWrVqpW8vb1ls9k0f/78Ii/H4XBo5syZatasmXx9fRUcHKyBAwfql19+KfmiK6CS2E6JiYkFvpeKs91x3u+//64ZM2aoR48eqlu3rux2u2rWrKnIyEh9+eWXRVqWu95LnCbBzdatW6eIiAj5+PhowIABqly5spYuXar+/fvr119/1ZgxYy65DIfDod69eyshIUFt27ZVZGSkkpKSNGfOHH322WfasmWLgoODS+HRWFNJbKNsvXv3VvPmzXNNDw0NLbmCK6Dnn39eycnJCgoKUq1atZScnFys5YwYMUJz5sxReHi4Ro0apT/++EOLFi3S6tWrtWXLFjVs2LCEK69YSmo7SVKnTp3UuXPnXNPzen+hcGbOnKkpU6aofv366tGjh4KDg5WUlKTly5dr+fLlWrBggfr371+oZbntvWTgNhkZGaZ+/frG29vb7Nixwzn9+PHjplGjRsZut5v9+/dfcjlz5841kszAgQONw+FwTp81a5aRZIYPH34lyq8QSmobzZs3z0gy8+bNu3LFVmBr1qxxbofY2NhiPddr1641kkzHjh3NuXPnnNNXrFhhJJkePXqUZMkVUklsp3Xr1hlJJjo6uuQLrOCWLl1qEhMTc03fsGGD8fLyMtWqVTNnz5695HLc+V5i95MbrV27Vnv37tX999/v8uuiatWqevbZZ5Wenq64uLhLLmf27NmSpNjYWJezl44YMUJhYWF699139ddff5V4/RVBSW0jXFndunVTSEjIZS0j+330j3/8Q3a73Tm9Z8+e6ty5s1avXq0DBw5c1joqupLYTrhy7r33XnXq1CnX9A4dOqhLly46duyYvv3220sux53vJUKNG2X3o+jRo0eu2yIiIiRJ69evL3AZZ8+e1Zdffqnrrrsu14eFzWZT9+7dlZaWpu3bt5dM0RVMSWyjnHbs2KFp06ZpypQpev/99/Xnn3+WSJ24fImJifL391e7du1y3VacbY0rKykpSTNmzFBsbKzi4+P1+++/u7skS/Py8pIkVap06V4r7nwv0afGjZKSkiQpz32LNWvWVEBAgHOe/Ozdu1cOhyPf/ZPZ05OSktShQ4fLrLjiKYltlNOrr77qct3X11fR0dEaN27c5RWKy5KWlqaDBw+qadOm8vT0zHV7zvcRyoYFCxZowYIFzuuVKlXSE088oZdffjnPbYjiO3DggD799FPVqlVLzZo1K3Bed7+XaKlxoxMnTkg6vysjL1WqVHHOcznLyDkfiqYktpEk1atXTzNnztTPP/+sM2fO6LffftN///tfVa9eXePHj9fMmTNLtG4UDe+j8iM4OFj/+te/9N133+n06dM6fPiwli9frgYNGmj69Ol65pln3F2ipWRkZGjQoEE6d+6cpkyZcsnA6O73EqEGKAWdOnXS448/roYNG8rX11d16tTRoEGDlJCQIB8fH8XExCgzM9PdZQJlXnh4uMaNG6fw8HD5+/urRo0a6t27t9atW6fg4GC9+uqrDGVRQhwOh4YMGaINGzYoKipKgwYNcndJl0SocaPsJJtfYj158mS+abcoy8g5H4qmJLZRQcLDw9W+fXulpqZq9+7dxV4OLg/vo/KvZs2a6t27tzIzM4s8pgpyczgcevjhh7VgwQI9+OCDeuONNwp1P3e/lwg1blTQvsVDhw7p9OnTlzyWPywsTB4eHvnunyyoTwgurSS20aUEBQVJOr8vGu7h7++vWrVqad++fcrKysp1O++j8oH3UslwOBwaOnSo4uLiNHDgQM2fP18eHoWLC+5+LxFq3Cj70LnVq1fnui0hIcFlnvz4+vrq5ptv1k8//ZRrICtjjNasWSN/f3+1atWqhKquWEpiGxUkKyvLeWQah7q6V6dOnZSWlqZNmzblui17W3fs2LG0y0IRZLfQMJhl8WUHmv/+97/q37+/4uPji9zx2q3vpSs2Ag4uKSMjw4SFhRU4sNu+ffuc0//44w+ze/duc/z4cZflMPjelVNS22j79u25lp2ZmWnGjh1rJJkuXbpcqYdQ4VxqULeUlBSze/duk5KS4jKdwfdKV3G3U17vJWOMmTFjhpFkGjZsaDIzM0u63AohKyvLDB482Egy9913n8nIyChw/rL4XiLUuNnatWuNl5eXqVy5somKijJPPfWUCQkJMZLM1KlTXebNfrFd/CGQlZVlIiIijCTTtm1bM27cOBMZGWlsNpupV6+eOXLkSCk+IuspiW0kydxwww3mwQcfNOPGjTNRUVGmUaNGRpK55pprzN69e0vxEVnP7NmzzeDBg83gwYNNixYtjCTTrl0757TZs2c7542Ojs53RNphw4YZSSY8PNw888wzZtCgQcZut5vq1aubn376qRQfkTWVxHYKCQkxDRo0MAMGDDBjx441I0eONDfddJORZAIDA82XX35Zyo/KOrKf84CAAPPcc8+Z6OjoXH85f9yVxfcSoaYM+PLLL83tt99uqlSpYnx9fc3NN99s3nvvvVzz5feFaYwxZ8+eNTExMaZ+/frGbrebmjVrmmHDhplDhw6VwiOwvsvdRmPGjDHt2rUzV199tfHy8jL+/v7mxhtvNM8//7xJTU0tpUdhXdnPe35/gwcPds5b0AdxVlaWeeWVV0x4eLjx9vY2V111lenfv7/Zs2dP6T0YCyuJ7fSvf/3LdOnSxdSuXdt4e3sbX19f07hxYzN69Gjz66+/lu4DsphLbZ+LP9vK4nvJZowxJbk7CwAAwB3oKAwAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACyBUAMAACzh/wGFp66NZGHWpgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x440 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "fit = posterior.sample(num_chains=4, num_samples=2000)\n",
    "y_pred_samples = fit['y_pred']\n",
    "\n",
    "# shape --> [num_draws, N_test]\n",
    "# Adjust y_test to match y_pred_samples shape (268, 8000) and compute RMSE per draw <-- Chatgpt did this because everything kep breaking, help\n",
    "rmse_samples = np.sqrt(np.mean((y_pred_samples - y_test[:, None])**2, axis=0))\n",
    "rmse_mean = rmse_samples.mean()\n",
    "rmse_ci = np.percentile(rmse_samples, [2.5, 97.5])\n",
    "\n",
    "idata = az.from_pystan(posterior=fit)\n",
    "print(az.summary(idata, var_names=['alpha', 'beta', 'sigma']))\n",
    "\n",
    "az.plot_forest(idata, var_names=['beta'], combined=True, colors=\"purple\")\n",
    "plt.axvline(0, linestyle='--', color='gray')\n",
    "plt.title(\"Posterior Distributions of Beta Coefficients\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}