experiments/eval_nobbs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a81aacda-6a13-44c8-ba56-d981566f0191",
   "metadata": {},
   "source": [
    "# Compare now-cast performance between our model and NobBS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a0776d32-7537-4542-811e-35ef8d5f9356",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import properscoring as ps\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set_style(\"whitegrid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "aeb4d6e4-1fce-4bf6-8eef-a905153ab785",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "plt.rc(\"text\", usetex=False)\n",
    "plt.rc(\"font\", family=\"serif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "006ae9e3-3cde-40a5-9cd6-7a9aa886e4d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_truth(ltla, report_date=\"2020-12-14\", n=14):\n",
    "    \"\"\"\n",
    "    Get the true count to compare the now-casts against.\n",
    "    \"\"\"\n",
    "    df = pd.read_csv(f\"../data/cases-{report_date}.csv\"\n",
    "                ).query(f\"`Area name` == '{ltla}' & `Area type` == 'ltla'\")\n",
    "    df = df.sort_values('Specimen date')\n",
    "    ys_sequence = df['Daily lab-confirmed cases'].values[-14:]\n",
    "\n",
    "    df = pd.read_csv(f'../data/cases-2020-12-21.csv')\n",
    "    df = df[(df['Area type']=='ltla')&(df['Specimen date']>='2020-10-09')&(df['Specimen date']<report_date)]\n",
    "    df['Specimen date'] = pd.to_datetime(df['Specimen date'])\n",
    "    df = df.sort_values(\"Specimen date\")\n",
    "    ltla_df = df[(df['Area name']==ltla)&(df['Specimen date']<report_date)]\n",
    "    truth = ltla_df[-n:]['Daily lab-confirmed cases'].values\n",
    "    return truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "90c78145-d81b-4d6a-b48a-d1d33d86f357",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/rjersakova/Documents/turing-rss/jbc-turing-rss-nowcasting\n"
     ]
    }
   ],
   "source": [
    "%cd .."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "cacc0fd4-d3b8-4409-bb48-db7246f01919",
   "metadata": {},
   "outputs": [],
   "source": [
    "from models.utils import empirical_credint\n",
    "from models.inference import Nowcaster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "914f332b-1053-4d93-ac10-aa1869b12eca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/rjersakova/Documents/turing-rss/jbc-turing-rss-nowcasting/experiments\n"
     ]
    }
   ],
   "source": [
    "%cd experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ad2a934d-b12e-4a80-92f5-1a9d1f9b3a33",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(f\"../data/cases-2020-12-14.csv\")\n",
    "all_ltlas = list(df.query(\"`Area type` == 'ltla'\")['Area name'].unique())\n",
    "all_ltlas.remove(\"King's Lynn and West Norfolk\")\n",
    "all_ltlas.remove(\"Kingston upon Hull, City of\")\n",
    "all_ltlas.remove(\"Eden\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a3f493b-21a7-4720-863b-9d7e93090a68",
   "metadata": {},
   "source": [
    "Compute absolute errors and CRPS for both models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5c8cca7c-c922-4d47-8709-f6c6c658d8cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "our_errors_all = []\n",
    "nobbs_errors_all = []\n",
    "\n",
    "our_crps_all = []\n",
    "nobbs_crps_all = []\n",
    "\n",
    "for ltla in all_ltlas:\n",
    "\n",
    "    truth = get_truth(ltla, n=18)\n",
    "\n",
    "    # original model\n",
    "    path = f\"../experiments/results/2020-12-14/{ltla.replace(' ', '_')}.pickle\"\n",
    "    model = pickle.load(open(path, 'rb'))\n",
    "    our_samples = np.array(model.x_smoothing_samples[-18:])\n",
    "    our_nowcast = our_samples.mean(axis=1)\n",
    "\n",
    "    # nobbs_nowcast\n",
    "    nobbs_df = pd.read_csv(f\"nobbs_results/{ltla}_2020-12-14_stan.csv\")\n",
    "    nobbs_samples = nobbs_df[[f\"sum_n_predict.{i}\" for i in range(1,19)]].T.to_numpy()\n",
    "    nobbs_nowcast = nobbs_samples.mean(axis=1)\n",
    "\n",
    "    # ABSOLUTE ERRORS\n",
    "    our_errors = np.abs(truth - our_nowcast)[-4:]\n",
    "    nobbs_errors = np.abs(truth - nobbs_nowcast)[-4:]\n",
    "    \n",
    "    our_errors_all.append(our_errors)\n",
    "    nobbs_errors_all.append(nobbs_errors)\n",
    "    \n",
    "    # CRPS\n",
    "    our_crps = [ps.crps_ensemble(truth[idx], x) for idx,x in enumerate(our_samples)][-4:]\n",
    "    nobbs_crps = [ps.crps_ensemble(truth[idx], x) for idx,x in enumerate(nobbs_samples)][-4:]\n",
    "    \n",
    "    our_crps_all.append(our_crps)\n",
    "    nobbs_crps_all.append(nobbs_crps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "599d4359-19ff-448b-adaa-97fd83fed1ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc(\"text\", usetex=False)\n",
    "plt.rc(\"font\", family=\"serif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "240243a9-5373-46e6-8ad0-0320092fe509",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAENCAYAAAASQ8TFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4kElEQVR4nO3de1jUZf7/8RcMAyikgqageAIVUXQ3ddNVt8zTfrMDlodCl7KviYdFXVzNw7Xmesi1lDy7nlZy2TQ1DDxruq6aFpZb39TQDPLH6noIRQ3kzPz+6HKuJlBnbIaZgefjuriuvD/33Pd7uIM37/l8PvfHw2QymQQAAAAA+Nk8nR0AAAAAAFQVFFgAAAAAYCcUWAAAAABgJxRYAAAAAGAnFFgAAAAAYCcUWAAAAABgJxRYQBU1YsQIdevWTeHh4c4OBQCA+zpy5IiioqIUHh6upUuX2vTaxYsXq0+fPgoPD1daWpqDIgSs4+XsAIDqLDMzU/Hx8bp06ZJu376tDRs2qH379hZ9/vjHP+rLL79Udna2mjRposmTJ6tr1673HXvNmjVaunSpli1bds9+6enpmjJlirKyslRcXKywsDBJ0vfffy9J6t27t8aNGyd/f3/za8rKypScnKzNmzeruLhYZWVlKi0tVcuWLdW9e3cNHDjQ1m8FAEDSmTNntGbNGp05c0ZeXl4qKSlRq1atNGLECLVp08bZ4TnUb37zG/3mN795oA8Gx48fry5duuill15yQGSAbTiDBThRaGioUlNT1bNnTxUXF2vChAnKzc216JOQkKA5c+YoMjJSqampVhVXtoiIiFBqaqoiIyNVv359paamKjU1Vf/85z/1+uuvKykpSW+88YbFa95++20tXrxY8+bNU0pKirZt26b169fr+vXrWrJkiV3jA4DqYseOHYqJiVGvXr20fft2paamavv27erVq5diYmL0wQcfODtEAFagwAJcRK9evfSf//xHM2bMcHYoZj169FDLli118OBBi/YtW7aoX79+5rNdklSvXj1NmTKlskMEgCrhzJkzmjJliqZMmaJ+/frJ0/OHP9E8PT319NNPa8qUKZo+fbpOnz7t5EgB3A8FFuAievfurZiYGO3YsUNbt269b/8DBw5owIAB6tu3r5544glNnTpV165dq7Dvl19+qZiYGPXp00e9e/dWSkqK1XGVlpaaE/0dJSUlunTpUrm+ERER2rx5s9VjAwB+sGLFCtWoUUPPPvtshcejoqLk6+urlStXSpLeeeedcvccHTt2TP369VN4eLhFHomKitKjjz6qnj176siRI4qJiTFfinfr1q1yc125ckVRUVF65JFHFBMTo61btyo6Olpdu3bVn/70JxUXF2v//v2KiYlR9+7d9frrr6uoqMhijBs3buj1119Xz5499dvf/lZRUVHatWtXubm2bt2q3r17q0+fPoqJidGXX35Z4fv/73//q/Hjx+uJJ54w58u79QWcjQILcCGvvfaa2rZtq9mzZ+vbb7+9a79du3Zp7NixGj16tPbt26c9e/bo6tWriomJ0e3bt8v1X79+vdasWaMPP/xQo0eP1uTJk3X06NF7xlJWVqYtW7YoIyNDr776qsWxzp07a9++fZo2bZrOnDljbvfw8FBQUJCN7xoAqreysjIdOXJEbdu2ldForLCPt7e32rRpo48++khlZWUaNmyY5syZY9Gna9euWr16dbnX3rkU/caNGzp8+LD+/ve/a9++fapXr16FczVo0MB86XhGRoZKSkq0ceNGJSUlKTk5WbNnz9bNmzeVlJSkxMREbdmyRampqebXFxUV6ZVXXlFmZqa2b9+uvXv3atKkSXrttdf0/vvvm/sdPnxYU6dO1ejRo/Xhhx9qzZo1euedd8rFk5OToyFDhqi4uFj79u3T/v371aNHD7300kvKyMiw5lsMVCoKLMCFeHt7a/HixfL09NSECRPKfSIoSSaTSW+99ZY6d+6s3r17S5J8fHw0adIkZWRkaNOmTeVe87vf/U6+vr6SpAEDBqhFixYVbn5x9epVRUVFKSoqSr/+9a81Z84cjR49Wi+//LJFvxkzZuhXv/qVkpOTFRUVpV69emnu3LkWxRYAwDo5OTm6ffu26tate89+9erV0+3bt5WTk/NA8+Tl5WnUqFHy8PBQjRo19N5771lsYFSRsrIyDRgwQJIUFhamFi1aaN++ferfv78kqWXLlgoLC9Mnn3xifk1qaqq++uorxcfHy8/PT5LUvXt39ejRQwsWLDDnthUrVig0NNQ8vq+vr4YOHVouhnfeeUeXLl3S5MmTzQXosGHDVKtWLa1Zs+aBvheAI1FgAS6mcePGmjNnjr766ivNnz+/3PHMzExdunRJ7dq1s2hv3bq1fHx8Kjwz1bJlS4t/t23bVidPnlRZWZlF+483uUhLS9PWrVt14MABDRkyRAUFBeZ+DRo00D/+8Q9t2bJFL7/8sgwGg9avX6+oqKhyn6gCAOzLw8PjgV5Xp04diyKucePG5S4B/6mQkBAZDAbzv2vXrl2urU6dOsrOzjb/+9ixY5JULk+1b99eOTk5Sk9PV2lpqU6ePFluZ8RWrVqVi+Hjjz9W3bp11bRpU3ObwWBQaGiojh8/fs/4AWdgm3bABT355JNKS0tTUlKSunXrpho1apiP3fnkslatWuVeV6tWrQo/2fzpJ5S1atVScXGxrl+/ftdLRKQfPq2cMGGCRo4cqS1btigmJsbiePv27dW+fXtNmzZNJ06c0MyZM5WUlKQnnnhC3bp1s+k9A0B1VadOHdWsWfOu99HekZ2dLT8/P9WpU+eB5rlzNskWP84/ksxnv37a9uMP7HJyclSjRg15e3tb9Ktdu7b5eE5OjkpKSsrlsoceeqhcDDk5OcrNzVVUVJRF+/fffy+TyWTzewIcjTNYgIuaOnWqwsPDNWXKFH333Xfm9oCAAEnSzZs3y73m1q1b5uM/9tOt32/evCmj0ajAwMD7xtG4cWNJ0jfffGNu27Ztm0pLSy36dezYUX/+858lSV999dV9xwUA/MBgMKh79+46ffq0iouLK+xTVFSkr776So899pj5rNOds0g/LjLy8/MdH/B9BAQEKD8/v9xl7nfyVkBAgAICAmQ0GsttslHRphsBAQF6+OGHzVdY/PhxIj/d5RZwBRRYgIvy8fHRokWLVFhYqL/85S/m9tDQUAUHB+vkyZMW/c+ePavCwsIKzxydO3fO4t+nT59Wu3bt7ntpiCRdvnxZkizOdC1atEjnz58v1/fOeBUVeQCAuxszZozy8/O1bdu2Co+npKSooKBAo0aNMrfd+ZDsx0XJvTZIqix3ntf4013+Tp48qYCAAEVERMhgMCgyMrLctvM/zVd3xrty5Ypu3Lhh0X7kyBEtWrTIrrED9kCBBbiw5s2ba9asWRbXtnt4eOi1115TWlqa9u/fL0kqLCzUggULFBoaqhdeeKHcOGvWrDHfQ5WcnKyMjAzFxcXdd/7r169rxYoVqlmzpvmG5jvmzp1rEdfVq1eVkJCg+vXrq2/fvg/ydgGg2oqIiNC8efM0b9487dq1y3zJXVlZmXbu3Kn58+frzTffVOvWrc2vadKkiYKCgsy5oKCgQNu3b3dK/D8WFRWlNm3aaOHChcrLy5P0w31ZBw8e1B//+EfzpYNxcXH69ttvlZycLOmH+CvatGLYsGGqX7++5s6daz4rdunSJb3xxhsW3w/AVXiYuHgVcJrMzEzFx8fr0qVLqlmzpho2bKgNGzaU6/enP/1J/+///T8lJSWZ2/bv368VK1YoNzdXRUVF6tKliyZNmmS+gXnEiBH66quvlJ2drcTERC1ZssRcEMXFxZkLpvT0dE2ZMkVZWVkqLi42Pzy4uLhYBQUFateunUaPHm2RxPbv36/9+/fr9OnTMhgMys/PV0lJiR599FH9/ve/V0hIiKO+ZQBQpaWnp2v16tU6e/asjEajiouL1bp1a8XGxlZYTKSlpWnOnDkqLCxUaGioXn75ZQ0bNkzBwcF67LHHNGvWLMXExOjs2bO6ffu2wsLCNHjw4Ap367ujoKBAL7zwgrKysiRJLVq00JYtWzRo0CDz5eIVtTVp0kSbNm2Sr6+vbty4oYSEBH300Ufy9vaWr6+vYmNj9dRTT1nMlZKSYt7Vtl69eho/fryGDRumevXqqU2bNuaC6/Lly5o/f74+++wzBQQEyMfHRy+//LL69esnSVq8eLF27NihrKwsNWnSRFFRUVZ9kAg4AgUWAAAAANgJlwgCAAAAgJ1QYAEAAACAnVBgAQAAAICdUGABAAAAgJ14OTsAR/viiy/k4+Pj7DB+lsLCQrd/D+6ONXANrINrqArrUFhYqF/+8pdW9SWPwF5YB9fAOriGqrAOd8slVb7A8vHxUUREhLPD+FnS09Pd/j24O9bANbAOrqEqrEN6errVfckjsBfWwTWwDq6hKqzD3XIJlwgCAAAAgJ1QYAEAAACAnVBgAQAAAICdUGABAAAAgJ1QYAEAAACAnVBgAQAAAICdUGABAAAAgJ1QYDlJZGSkPDw8rPpq06aN1X0jIyOd/dYAAACAaosCy0lOnTolk8lk1VfTyTus7nvq1ClnvzUAAACg2qLAAgAAAAA7ocACAAAAADuhwAIAAAAAO6HAAgAAAAA7ocACAAAAADuhwAIAAAAAO6HAAgAAAAA7ocACAAAAADvxcubke/fu1d///nfl5+fr5s2bql27tl566SX1799fklRcXKzly5drz5498vLykr+/vyZOnKhOnTo5M2wAAAAAqJDTzmC98847WrlypRISErR161bt2bNHzZs31yeffGLuM3v2bO3evVsbNmzQjh07NGDAAA0fPlzp6enOChsAAAAA7sopBdaFCxe0YMECzZw5U0FBQZIko9Go1157TUOHDpUkZWZmavPmzRoxYoQCAwMlSYMGDVJISIgWLlzojLABAAAA4J6cUmClpqaqVq1aat++vUV7gwYN1K5dO0nS/v37ZTKZ1LlzZ4s+Xbp00dGjR5WXl1dp8QIAAACANZxyD9bnn3+uRo0aae/evVq/fr2uX7+uOnXqaODAgRo4cKAk6ezZs/L09FTDhg0tXhsSEqKSkhJlZGSUK9AqUlhYWCUuKawK78GdFRQUsAYugHVwDdVtHapCHqlua+aqWAfXwDq4hqq8Dk4psC5fvqyLFy9q3bp1WrZsmerWrau9e/fqj3/8o7777juNHj1aOTk58vX1lcFgsHitv7+/JCknJ8equXx8fBQREWH391C5MqvAe3Bv6enprIELYB1cQ1VYB1uSelXII1VhzaoC1sE1sA6uoSqsw91yiVMuESwsLNTt27c1adIkPfzww/L09NSTTz6pXr16adWqVcrPz3dGWAAAAADwszilwPLz85OkclVrRESE8vPz9c033yggIEAFBQUqLS216JObmytJCggIqJxgAQAAAMBKTimwQkNDJUkmk8mi/c7lgCaTSeHh4SorK9OlS5cs+ly4cEFeXl4KCwurnGABAAAAwEpOKbCeeOIJSdKZM2cs2r/++mv5+vqqRYsW6t27tzw8PHT8+HGLPmlpaerWrZv5LBgAAAAAuAqnFFj9+vVTu3bttHjxYvN265999pn27t2rUaNGqWbNmgoNDdXgwYO1atUqXb9+XZKUnJysrKwsxcfHOyNsAAAAALgnp+wiaDAYtHbtWi1YsEBPP/20vL295e3trddff12DBw8295s+fbqWLVum6OhoGY1G+fn5ad26dW6/4wgAAACAqskpBZYk1alTR3PmzLlnH6PRqPj4eM5YAQAAAHALTrlEEAAAAACqIgosAAAAALATCiwAAAAAsBMKLAAAAACwEwosAAAAALATCiwAAAAAsBOrC6y4uDjFxcXp8uXLjowHAAAAANyW1QXW4cOH1b9/f9WrV8+R8QAAAACA27K6wGrdurV69+4tL6+Kn0185coVuwUFAAAAAO7I6gKrS5cu+vTTT+96fNSoUXYJCAAAAADcVcWnoypgMBg0adIktW7dWqGhoapZs6bF8e+++87uwQEAAACAO7G6wPrrX/8qSbp8+bL+9a9/lTvu4eFht6AAAAAAwB1ZXWC1bt1aKSkpdz3ev39/O4QDAAAAAO7L6gJr3Lhx9zz+pz/9yaaJL1y4oGeeeUZNmjQpdywpKUm1atWSJBUXF2v58uXas2ePvLy85O/vr4kTJ6pTp042zQcAAAAAjmZ1gdWzZ0/zf1+/fl05OTkKCAhQYGCgJD1QwRMZGamkpKR79pk9e7bS0tK0ceNGBQYGasuWLRo+fLjee+89RURE2DwnAAAAADiK1bsIStKxY8cUFRWlbt266emnn1a3bt3Uv39/ffzxxw4JLjMzU5s3b9aIESPMhdygQYMUEhKihQsXOmROAAAAAHhQVhdYH3/8sUaMGCGDwaDBgwcrNjZWgwcPlqenp0aMGKFPPvnE7sHt379fJpNJnTt3tmjv0qWLjh49qry8PLvPCQAAAAAPyupLBJcsWaLZs2fr+eefL3csJSVFixcvVpcuXWyaPDs7WxMnTtT58+f1/fffKzIyUrGxsQoPD5cknT17Vp6enmrYsKHF60JCQlRSUqKMjAy1b9/epjkBAAAAwFGsLrCuXbtWYXEl/bCD4IoVK2ya2GAwyMvLS8OGDVNkZKRyc3M1e/ZsDR48WElJSWrfvr1ycnLk6+srg8Fg8Vp/f39JUk5Ozn3nKSwsVHp6uk2xuaKq8B7cWUFBAWvgAlgH11Dd1qEq5JHqtmauinVwDayDa6jK62B1gVVSUiKTyVTh865KS0tVXFxs08TBwcHavn27+d/+/v6aOXOmDh06pIULFyoxMdGm8e7Gx8enCmyGkVkF3oN7S09PZw1cAOvgGqrCOtiS1KtCHqkKa1YVsA6ugXVwDVVhHe6WS6y+B+uRRx7RxIkTdeXKFYv2y5cva9KkSerQocPPi1CSr6+vWrVqpS+++EKSFBAQoIKCApWWllr0y83NNR8HAAAAAFdh9RmsSZMm6cUXX9QTTzyhevXqqXbt2rpx44auXbumoKAgbdiwwaaJv//+e/n4+Mjb29ui3WAwqKysTJIUHh6uHTt26NKlSwoJCTH3uXDhgry8vBQWFmbTnAAAVBWRkZE6ffq03cdt27atTp06ZfdxAaC6sPoMVlBQkFJSUjRixAgFBwersLBQDRs2VGxsrLZu3aqgoCCbJn7jjTe0b98+i7aioiJ9/fXXatOmjSSpd+/e8vDw0PHjxy36paWlqVu3bvLz87NpTgAAqopTp07JZDJZ9dV08g6r+1JcAcDPY/UZrAMHDkiSRo4cqfj4eLtM/re//U2PPvqo6tevr9LSUr311lu6fv263nrrLUlSaGioBg8erFWrVqlHjx4KDAxUcnKysrKyNH/+fLvEAAAAAAD2YnWB9fvf/169evXSI488opo1a/7siV955RVt2rRJr776qqQfdgQMCwtTYmKixXbv06dP17JlyxQdHS2j0Sg/Pz+tW7fO7W+KAwAAAFD1WF1ghYaGavny5XabODw8XK+//vp9+xmNRsXHx9vtrBkAAAAAOIrV92A1adLEvHtfRaZNm2aXgAAAAADAXVl9BmvYsGGaOHGinn/+eYWGhpa7TPDLL7+0e3AAAAAA4E5sKrAk6dChQ46KBQAAAADcmtUFVnBwsMaNG1fhMZPJpGXLltktKAAAAABwR1YXWN26ddNzzz131+PffPONXQICAAAAAHdl9SYXNWvW1F/+8hddu3atwuOvvfaa3YICAAAAAHdkdYH1j3/8Q8HBwfLz83NkPAAAAADgtqy+RLB169bmjS4qYjKZ5OHhYY+YAAAAAMAtWX0GKzIyUhkZGXc9/vzzz9slIAAAAABwVzadwRo7dqy6du2q5s2bl7tU8ObNm3YPDgAAAADcidUF1qxZsyRJmZmZFR7n8kAAAAAA1Z3VBVZYWJhWr15d4TGTyaSRI0faLSgAAAAAcEdWF1gxMTFq1KjRXY///ve/t0tAAAAAAOCurN7k4sUXX7zn8X79+v2sQIYMGaLw8HBduHDhZ40DAAAAAM5i9RksSfr222+1evVqHT9+XJJ04MABLV68WBEREerbt+8DB7F3716dOHGiwmN5eXlKSEjQ0aNHZTAYFBQUpKlTp6ply5YPPB8AAAAAOILVZ7C+/PJLPf/88zp27JiaNGlibu/QoYMWLlyovXv3PlAARUVFSkhI0OOPP17h8fHjxys9PV0pKSnatWuX2rdvr5iYGF25cuWB5gMAAAAAR7G6wEpISNDYsWN18OBBJSYmqlatWpKk3/zmN1q3bp0SExMfKIANGzaoXbt2ateuXbljR48e1ZEjRzRu3DjVqFFDkjRmzBiVlZVp5cqVDzQfAAAAADiK1QXWpUuX9L//+7/y9Cz/kuDgYBUVFdk8+Y0bN7R27VpNmDChwuN79+6V0WhUx44dzW3e3t7q0KGD9u3bZ/N8AAAAAOBIVhdYJSUlKisrq/BYcXGxcnJybJ58+fLlevbZZ++6O+HXX3+t+vXry9vb26I9JCRE2dnZunbtms1zAgAAAICjWL3JRfv27TVu3DhNnjxZjRs3NrffunVLc+fOVYcOHWya+Pz589q9e7d279591z45OTny8/Mr136n7caNG6pbt+495yksLFR6erpNsbmiqvAe3FlBQQFr4AJYB9dQ3daBPAJ7qW4/O66KdXANVXkdrC6wJk+erOjoaPXt21d169ZVbm6u+vTpo8uXL6t+/frasGGDTRMvWLBAsbGxeuihh2wO2hY+Pj6KiIhw6ByOl1kF3oN7S09PZw1cAOvgGqrCOtiS1MkjsJeq8LNTFbAOrqEqrMPdconVlwgGBwcrJSVFsbGxatSokerXr6/AwEC9+uqr2rp1qxo0aGB1MJ999pnOnTun6Ojoe/YLCAhQXl5eufY7bXXq1LF6TgAAAABwNJueg1WnTh3Fx8f/7EmPHj2q0tJSDRw40NyWnZ0tSYqNjZXRaNSECRPUqlUrnTp1SkVFRRb3YV24cEH16tW77+WBAAAAAFCZbCqw7GX8+PEaP368RdvSpUu1bNkyrV69WiEhIZIkg8GgTZs26fPPP1fnzp0l/fDcrH//+9966qmnKj1uAAAAALgXqy8RdIbu3bure/fuWrx4sfLz8yVJK1eulKenp0aNGuXk6AAAAADAklPOYP3YoUOH9Pbbb5e7RDA1NVWStGTJEi1YsED9+/eXp6engoKClJSUZNM9XwAAAABQGZxeYD3++ON6/PHH73rcz89PM2bMqMSIAAAAAODBuPQlggAAAADgTmwusAoKCvTpp5/qwIEDkn54GDAAAAAAwMZLBFesWKG1a9cqPz9f9erVU69evTRjxgyVlJTo7bfflq+vr6PiBAAAAACXZ/UZrMTERCUlJWno0KGaN2+eatWqJUlasGCBGjVqpMWLFzssSAAAAABwB1afwdqyZYveffddhYaGSpLeeecdSZK3t7cmT56sQYMGOSRAAAAAAHAXNt2Ddae4+ikvLy8VFxfbJSAAAAAAcFdWF1ilpaX69ttvKzx2/vx5CiwAAAAA1Z7Vlwg+99xzio6O1pAhQ/TII4+ooKBAJ06c0JkzZ/S3v/1NQ4cOdWScAAAAAODyrC6wRo4cqcuXL+uvf/2rJMlkMul3v/udJGno0KEaPny4YyIEAAAAADdhdYHl4eGhP//5z3rllVf08ccfKycnRwEBAeratauaNGniyBgBAAAAwC1YXWBlZGQoLCxMTZs2VdOmTc3tZ86c0YoVKxQdHa2AgACHBAkAAAAA7sDqTS4mTpxYYbvRaFRmZqbi4+PtFhQAAAAAuCOrCyyTyVRhe1hYmBYsWKCcnBy7BQUAAAAA7uielwh++umnSktLkyRlZ2dr+fLlFRZaly9fVn5+vmMiBAAAAAA3cc8CKy0tTcuWLZP0wyYXS5cuLT+Al5caNWqkyZMnWz1pVlaWNm7caC7e8vLyVLduXcXGxqpHjx7mfsXFxVq+fLn27NkjLy8v+fv7a+LEierUqZPVcwEAAABAZblngRUXF6e4uDhJUv/+/ZWSkmKXSQ8fPqydO3cqKSlJTZs2VVlZmRYsWKDRo0dr/fr1evTRRyVJs2fPVlpamjZu3KjAwEBt2bJFw4cP13vvvaeIiAi7xAIAAAAA9mL1PVizZs2y26QNGjTQ2LFjzbsRenp6atSoUSorK9OBAwckSZmZmdq8ebNGjBihwMBASdKgQYMUEhKihQsX2i0WAAAAALAXqwus9u3b3/P4ihUrrJ60T58+GjRokEVbbm6uJJmLqf3798tkMqlz584W/bp06aKjR48qLy/P6vkAAAAAoDJY/RysTz/99J7HU1NTNWbMmAcK4sqVK5o1a5batm2roUOHSpLOnj0rT09PNWzY0KJvSEiISkpKlJGRcd+iT5IKCwuVnp7+QHG5kqrwHtxZQUEBa+ACWAfXUN3WgTwCe6luPzuuinVwDVV5HawusGJiYuTh4WHXybOysjR8+HBlZWXpscce0/Lly+Xv7y9JysnJka+vrwwGg8VrfnzcGj4+PlXgfq3MKvAe3Ft6ejpr4AJYB9dQFdbBlqROHoG9VIWfnaqAdXANVWEd7pZLrC6wmjRpojlz5li05efnKzMzUwcOHFBMTIzNQTVp0kQffvihcnNztXjxYkVFRWnFihXsEggAAADALVldYI0ZM8a8u9+PPf7444qKitIbb7yh3/72tw8UhL+/v6ZNm6aPPvpIM2fO1Pbt2xUQEKCCggKVlpZanMW6c69WQEDAA80FAAAAAI5i9SYX/fv3v+uxwMBAnT171upJCwoKyj2w2MPDQ61atdK5c+dUVFSk8PBwlZWV6dKlSxb9Lly4IC8vL4WFhVk9HwAAAABUBqsLrP/+97/lvi5evKj09HQtWbJExcXFVk86YsQIffHFF+XaL168KH9/fxmNRvXu3VseHh46fvy4RZ+0tDR169ZNfn5+Vs8HAAAAAJXB6ksEe/bseddNLnx9fTVv3jybJl66dKkSEhIUEBAgk8mkpKQknTx5UmPGjJGHh4dCQ0M1ePBgrVq1Sj169FBgYKCSk5OVlZWl+fPn2zQXAAAAAFQGqwus4OBgjRs3zqLN09NTdevWVbt27VS7dm2rJ/3DH/6g999/XzExMTIYDCosLFSdOnU0f/58PfPMM+Z+06dP17JlyxQdHS2j0Sg/Pz+tW7fO7XccAQAAAFA1WV1gPfnkk3ruuefsMmnHjh3VsWPH+/YzGo2Kj49XfHy8XeYFAAAAAEey+h6s11577Z7HDx8+/LODAQAAAAB3ZnWBdT8LFy6011AAAAAA4Jbueokg9zkBAAAAgG3uWmDVrVtXL774olWDmEwmbdq0yW5BAQAAAIA7umuB1aZNG8XFxVk90MmTJ+0SEAAAgLuIjIzU6dOn7T5u27ZtderUKbuPC8Dx7noP1urVq20ayNb+AAAA7u7UqVMymUxWfTWdvMPqvhRXgPuyept2SSooKNB7772nY8eOKScnRwEBAerWrZteeOEF+fr6OipGAAAAAHALVhdY169fV0xMjDIyMuTj46NatWrp66+/1uHDh7Vlyxb9/e9/V2BgoCNjdWm/mLlPN/OLHTZ+syk77T5m7RpG/d+MvnYfFwAAAKiurC6w5s+fr+DgYL399tsKDw83t589e1YJCQlasGCB5s6d65Ag3cHN/GKdn/eUQ8ZOT093yK6OjijaAAAAgOrM6gIrLS1Ne/bskbe3t0V7eHi4li5dqieffNLuwQEAAACAO7H6QcNGo7FccXWHj4+PjEaj3YICAAAAAHdkdYH10EMP6eDBgxUeO3DggPz9/e0WFAAAAAC4I6svERw9erTi4uLUpUsXRUZGqnbt2rp586ZOnjyptLQ0LV261JFxAnbHs0sAAABgb1YXWL169dL8+fM1f/58HT161NweHBys+fPnq2fPnlZPmp6ernfffVcnTpyQl5eXSktL1bVrV40ZM8ZiJ8K8vDwlJCTo6NGjMhgMCgoK0tSpU9WyZUur5wLuxpYiqNmUnQ7bxAQAAABVh03PwerXr5/69eunzMxM83OwQkNDbZ40Pj5eLVu2VHJysmrWrKkrV67o5Zdf1pEjR5Sammp+ptb48eOVl5enlJQU1ahRQ4sWLVJMTIxSU1PVoEEDm+cFAAAAAEey+h6sHwsNDVXHjh0VGhqqW7duPdDEEydOVM2aNSVJDRo00PDhw3X+/HkdOnRIknT06FEdOXJE48aNU40aNSRJY8aMUVlZmVauXPlAcwIAAACAI1ldYB0/flxTp07Vzp0/PDspIyNDPXv2VOfOnfX888/r8uXLVk+6bds2NW3a1KKtfv36kmQu2Pbu3Suj0aiOHTua+3h7e6tDhw7at2+f1XMBAAAAQGWxusBKSkpSUVGRWrduLUmaPXu28vLyNHnyZDVt2lTz58+3etKKtns/f/68PDw81KlTJ0nS119/rfr165frGxISouzsbF27ds3q+QAAAACgMlh9D1ZWVpa2bt0qg8Ggy5cv65NPPtGMGTMUHR2tkpIS9e/f/4GDKC0t1fvvv6+BAweqefPmkqScnBz5+fmV63un7caNG6pbt+59xy4sLFR6evoDx2YLR81TUFDgsLEr63tTFfC9cj5H/izAetVtHSozjzhSVXgPVQHr4HzV7XeYq6rK62B1gWUwGGQwGCRJe/bska+vr5555pkfBvHy+lkPGl6xYoWMRqOmTZv2wGPcjY+PjyIiIuw+bnmZDpsnPT3dQWM7Luaqh++VK3DczwJsURXWwZakXnl5xJH4HeYaWAdXUBV+h1UFVWEd7pZLrL5E0MvLS59//rm+++47JSUlqXfv3uaHC1+7dk0lJSUPFFhycrJ2796tNWvWmDe9kKSAgADl5eWV63+nrU6dOg80HwAAAAA4itVnsMaMGaOXXnpJJSUl8vHx0ciRIyVJW7Zs0dq1a/XrX//a5slTUlKUmJio9evXl7vcr1WrVjp16pSKioos7sO6cOGC6tWrZ9XlgQAAAABQmaw+g9WjRw/t2rVLixYt0q5du8wP+w0JCdGoUaP06quv2jRxamqq1q5dq8TERD388MOSpIMHD2rTpk2SpL59+6q4uFiff/65+TVFRUX697//rb59+9o0FwAAAABUBpseNNy4cWM1btzYou1Bzlxt27ZN06dP1/jx43Xs2DFz+4kTJ8zFVvfu3dW9e3ctXrxYf/vb31SjRg2tXLlSnp6eGjVqlM1zAgAAAICj2VRgZWRkaNWqVfriiy909epV1a9fX4888ohGjhyp0NBQq8eZM2eOCgsL9dZbb5U7FhcXZ/7vJUuWaMGCBerfv788PT0VFBSkpKQkNWjQwJawAQAAAKBSWF1gHTlyRKNHj5afn59atmypFi1a6ObNmzp48KB27dqllStXqlu3blaNdfz4cav6+fn5acaMGdaGCAAAAABOZXWB9dZbb2nUqFGKjY212HSiqKhIq1at0rx587R9+3aHBAkAAAAA7sDqTS7y8vIUFxdnUVxJkre3t8aOHavc3Fy7BwcAAAAA7sTqM1hBQUEymUzy8PAod6ysrIxt0wEA+Jl+MXOfbuYXO2TsZlN22n3M2jWM+r8Z7OwLAD9mdYE1cuRI/fnPf9a4ceMsiqnvvvtOixYt0iuvvOKQAAEAqC5u5hfr/Lyn7D5uenq6IiIi7D6uI4o2AHB3dy2wevXqVa4tJydHmzdvVq1ateTv76/c3FzdunVLRqNRaWlpeuop+ycFAAAAAHAXdy2wcnNz1bNnT6sGMZlMOnjwoN2CAgAAAAB3dNcCKzg4WH/5y1+sHuiZZ56xS0AAAAAA4K7uWmClpKRYNUBZWZkOHTqknJwce8UEAAAAAG7J6k0ufurMmTNKSUnRjh07lJ2dXeHuggAAAABQndhUYF2/fl3btm1TSkqKzp49K4PBoF/96ld69dVXtWHDBkfFCAAAAABu4b4FVnFxsf75z3/qgw8+0EcffaTS0lJ17NhRwcHB2r59u/z8/CRJBoPB4cECAAAAgCu7Z4E1c+ZM7dq1Szdv3lStWrUUHR2tF198UWFhYXruuefMxZUkxcTEODxYAAAAAHBl9yyw3nvvPRkMBk2ePFlDhw6Vt7d3ZcUFAAAAAG7H814H//Wvf2ncuHHatm2bpk6dqk8//dSuk1+9elXDhw9XeHi4XccFAAAAAGe45xmsBg0aKDY2VrGxsTp58qRSUlK0aNEi9e3bVyUlJRZ9z58/r2bNmlk98b59+zRv3jx5ed37NrDi4mItX75ce/bskZeXl/z9/TVx4kR16tTJ6rkAAAAAoDLc8wzWj7Vr107Tp0/XO++8o4YNG6pZs2aKi4vTe++9pytXrig+Pt6midesWaN169apQ4cO9+w3e/Zs7d69Wxs2bNCOHTs0YMAADR8+XOnp6TbNBwAAAACOZnWBdYfRaFSfPn20dOlSzZo1S4WFhRo9erTOnTtn0zgbN2687xmvzMxMbd68WSNGjFBgYKAkadCgQQoJCdHChQttDR0AAAAAHMrmAuvHAgMD9fLLL2vr1q0KCwuz6bX3uzRQkvbv3y+TyaTOnTtbtHfp0kVHjx5VXl6eTXMCAAAAgCP9rALrx9avX2+voczOnj0rT09PNWzY0KI9JCREJSUlysjIsPucAAAAAPCg7n8ayUp16tSx11BmOTk58vX1LfcQY39/f/Px+yksLKy0+7UcNU9BQYHDxq5q97IN2nheuUVlDhm72ZSddh/T39tTW6Kb2X3cqsqRPwuwXnVbh8rMI5Jjfi+TR1wH3y/HePbZZ/XNN9/YfdwWLVpo27Ztdh8XVTuX2K3AclU+Pj6KiIhw+Dw1m4/V4ONXHD6PPdVs3kAREfudHYZd5RZl6vy8p+w+bnp6ukP+P2o2ZWel/P9ZVThqHWCbqrAOtiT1ysojP8h0yFyOWzPHxFt18f1yFFv2Amg2ZadD/laAbapyLnHpAisgIEAFBQUqLS21OIuVm5trPu4qbn8b77AfVkf+cQ84W2RkpE6fPm33cdu2batTp07ZfVwAAIB7sds9WI4QHh6usrIyXbp0yaL9woUL8vLysnljDQCu59SpUzKZTFZ9NZ28w+q+FFcAAMAZXPoMVu/evfX222/r+PHjCgkJMbenpaWpW7du8vPzc2J0AACgKvrFzH26mV/skLEdcfVI7RpG/d+MvnYfF8CDcekCKzQ0VIMHD9aqVavUo0cPBQYGKjk5WVlZWZo/f76zwwMAwK5qNl+oduunOGbw4/YfsmbzBpKq3r0sN/OL3e5+XgCuw2kF1ptvvqljx46ZL/+LioqSJG3ZskXe3t7mftOnT9eyZcsUHR0to9EoPz8/rVu3zu1vigMA4KccdT8vf9gDQOVxWoE1efJkq/oZjUbFx8crPj7ewREBAAAAwM/j0ptcAAAAAIA7ocACAAAAADuhwAIAAADws0VGRsrDw8OqrzZt2ljdNzIy0tlvzSYUWAAAAAB+Np5t+QMKLAAAAACwEwosAAAAALATCiwAANfNAwBgJ057DhbgCDWbL1S79VMcM/hx+w9Zs3kDSfZ/qChgK1uub282ZadDHoYLuApyCYCfgwLLjhz7RPtMu49Yu4bR7mM62+1v4x3yh196eroiIiLsPq5j/58B4I4c93uBPGItcolr+MXMfbqZX+yQsR3xPatdw6j/m9HX7uPC/VBg2YkjP83l02K4G5Ii8GAc9buePAJ3dDO/mEIXbokCC4DdkRQBAEB1xSYXAAAAAGAnFFgAAAAAYCdcIogqhxvEAQA/F7kEwINy+QLr2rVrmjt3rnkL4fDwcE2bNk1BQUFOjgyuiBvEXQNbHLsGNhsBHgy5xDWQS1wDucR2Ll1gFRUV6ZVXXlGzZs20c+dOeXh4aNq0aXrppZf0wQcfyM/Pz9khAqgAWxy7BjYbAeDOyCWugVxiO5cusFJSUnT27FktX75cXl4/hDpx4kQ99thj2rhxo1599VUnRwjgbri8xvn49BeAuyOXOB+5xHYuXWDt3btXDRs2VOPGjc1tDz/8sFq0aKF9+/a5dYEVGRmp06dPW93f403r+rVt29Z8OSXujTVwHC6vcQ23v413dgg2qap/nDgKv8NcA+vgOOQS10AusZ2HyWQyOTuIu+nevbuaN2+upKQki/ZRo0bpk08+0RdffHHfMb744gv5+Pg4KMLKUVBQIF9fX2eHUa2xBo7z7LPP6ptvvrH7uC1atNC2bdvsPm5VVR3XwdpLU8gjsBfWwXGq4+8wV1Qd16GiXOLSZ7Bu3LhR4X1W/v7+ys/Pt+oXlY+Pj0Ou76xMjrpGFdZjDRzn3LlzVvdlHRynuq1Denq61X3JI7AX1sFxqtvvMFdV3dbhbrmE52ABAAAAgJ24dIFVp04d5eXllWvPzc1VjRo1OM0OAAAAwKW4dIEVHh6uixcvlmu/ePGiWrVq5YSIAAAAAODuXLrA6tu3ry5evKgLFy6Y27Kzs5WRkaG+fXkYJQAAAADX4tIF1nPPPadWrVppwYIFKikpUVlZmRISEhQSEqLo6GhnhwcAAAAAFly6wPL29lZiYqIMBoOeeuopPfnkk8rNzdX69esr3F0QAAAAAJzJpbdpl6R69eopISHB2WEAAAAAwH259IOG7aEqPCASAGBfhYWF+uUvf2lVX/IIAKAid8slVb7AAgAAAIDK4tL3YAEAAACAO6HAAgAAAAA7ocACAAAAADuhwAIAAAAAO6HAAgAAAAA7ocACAAAAADuhwHJRV69e1fDhwxUeHu7sUACgnCFDhig8PFwXLlxwdii4B3IJAFdVlfOIl7MDQHn79u3TvHnz5OXF8jhLenq63n33XZ04cUJeXl4qLS1V165dNWbMGAUGBjo7vGojKytLGzduVFpamiQpLy9PdevWVWxsrHr06OHc4KqxvXv36sSJE84OA/dBLnE+colrIJe4nqqeRziD5YLWrFmjdevWqUOHDs4OpdqKj4/XzZs3lZycrO3btysxMVEfffSRoqOjVVBQ4Ozwqo3Dhw9r586dWrhwobZu3ardu3frl7/8pUaPHq3jx487O7xqqaioSAkJCXr88cedHQrug1zifOQS10AucS3VIY9QYLmgjRs3qlmzZs4Oo9qbOHGiatasKUlq0KCBhg8frvPnz+vQoUNOjqz6aNCggcaOHaumTZtKkjw9PTVq1CiVlZXpwIEDTo6uetqwYYPatWundu3aOTsU3Ae5xDWQS5yPXOJaqkMe4boBF8TlHM63bds2eXt7W7TVr19fknTr1i1nhFQt9enTp1xbbm6uJHF5jRPcuHFDa9eu1aZNm7R161Znh4P7IJc4H7nENZBLXEd1ySOcwQIq8NOEKEnnz5+Xh4eHOnXq5ISIIElXrlzRrFmz1LZtWw0dOtTZ4VQ7y5cv17PPPqtGjRo5OxTALZBLXBO5xHmqSx6hwAKsUFpaqvfff18DBw5U8+bNnR1OtZOVlaU+ffroscceU2lpqZYvXy5/f39nh1WtnD9/Xrt379bo0aOdHQrgtsglzkUuca7qlEcosAArrFixQkajUdOmTXN2KNVSkyZN9OGHH+rEiRNq1qyZoqKi9Nlnnzk7rGplwYIFio2N1UMPPeTsUAC3RS5xLnKJc1WnPEKBBdxHcnKydu/erTVr1phvVIZz+Pv7a9q0aapbt65mzpzp7HCqjc8++0znzp1TdHS0s0MB3Ba5xHWQSypfdcsj3AEL3ENKSooSExO1fv161a1b19nhVDsFBQXy8fGRh4eHuc3Dw0OtWrXS3r17VVRUVOE9DrCvo0ePqrS0VAMHDjS3ZWdnS5JiY2NlNBo1YcKEKr3lLvBzkEuci1zifNUtj1BgAXeRmpqqtWvXKjExUQ8//LAk6eDBg7p69apeeOEFJ0dXPYwYMUITJkzQI488YtF+8eJF+fv7y2g0Oimy6mX8+PEaP368RdvSpUu1bNkyrV69WiEhIU6KDHB95BLnI5c4X3XLIxRYQAW2bdum6dOna/z48Tp27Ji5/cSJE+YEicqxdOlSJSQkKCAgQCaTSUlJSTp58qTGjBlj8WkkALgaconrIJegMnmYTCaTs4OApTfffFPHjh3TpUuXdPPmTbVu3VqStGXLFk5hV5JHH31UN2/erPBYXFycxo4dW8kRVU8nTpzQ+++/r5MnT8pgMKiwsFB16tTRkCFD9Mwzz5AUneDQoUN6++23lZ2drezsbIWFhcloNCo1NdXZoeEnyCXORy5xDeQS11Id8ggFFgAAAADYCbsIAgAAAICdUGABAAAAgJ1QYAEAAACAnVBgAQAAAICdUGABAAAAgJ1QYAEAAACAnVBgAQAAAICdUGABAAAAgJ1QYAFu5tq1a4qKitKjjz6q8PBwRUVFKTk52dlhAQDcBHkEcCwPk8lkcnYQAGw3ZcoUffDBBzp79qyzQwEAuCHyCOAYnMECAAAAADvxcnYAABxn5cqV2r9/v8rKylRcXKyWLVtq8uTJatCggUW/rVu3asWKFfLw8FC9evU0btw4DRs2TPXq1VObNm20Zs0aJ70DAIAzkUcA21FgAVXY2rVr9Y9//EOtW7dWaWmp3nzzTY0aNUrvv/++DAaDJOnw4cOaOnWq5s6dqwEDBqigoEB/+MMfJEkvvviixo4d68R3AABwJvIIYDsuEQSqsE2bNql169aSJIPBoBdeeEFfffWVTp06Ze6zYsUKhYaGasCAAZIkX19fjRgxwinxAgBcC3kEsB1nsIAq7NatWxozZoyysrJkMBhUXFwsSfrPf/6jX/ziFyotLdXJkyf1P//zPxava9WqlTPCBQC4GPIIYDvOYAFV1JkzZ/TSSy+pSZMmSklJUWpqqlavXi1JKioqkiTl5OSopKREtWrVsnjtQw89VOnxAgBcC3kEeDAUWEAVVFRUpF27dqmoqEijRo2Sl1fFJ6sDAgJkNBp169Yti/af/hsAUL2QR4AHxyWCQBU0fPhw5eXlSZI8PDzM7dnZ2Rb9DAaDIiMjdfr0aYv2c+fOOT5IAIDLIo8AD44zWEAVNWnSJElSYmKipB8+jaxom9y4uDh9++23Sk5OliQVFBTo3XffrbxAAQAuiTwCPBgPk8lkcnYQAKz33Xff6YUXXlBOTo5u376tRo0aVdhn7dq1On/+vNauXStPT08FBQXpscce01tvvaXg4GA9/fTTmjhxoiQpJSVFy5Ytk4eHh4KCgjRp0iQNGjRIY8eOVVxcXGW/RQCAA5FHAMeiwAJQzo0bN9S5c2fNmDFDQ4YMcXY4AAA3Qx5BdcYlgkA199lnn2njxo0WbXeunW/ZsqUzQgIAuBHyCGCJAguo5m7cuKE1a9bo6tWrkqTc3FwtWbJEnTp10q9+9SsnRwcAcHXkEcASuwgC1Vx4eLg6dOigmJgYeXt769atW+rSpYsmT57s7NAAAG6APAJY4h4sAAAAALATLhEEAAAAADuhwAIAAAAAO6HAAgAAAAA7ocACAAAAADuhwAIAAAAAO/n/dSC0LK9BzQIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(12, 4), sharey=True)\n",
    "\n",
    "pd.DataFrame(np.flip(np.array(our_errors_all))).boxplot(showfliers=False, whis=(5,95), ax=ax[1])\n",
    "pd.DataFrame(np.flip(np.array(nobbs_errors_all))).boxplot(showfliers=False, whis=(5,95), ax=ax[0])\n",
    "\n",
    "for i in range(2):\n",
    "\n",
    "    ax[i].set_xticklabels([1,2,3, 4])\n",
    "    ax[i].set_xlabel(\"Lag\", fontsize=16)\n",
    "    plt.setp(ax[i].get_xticklabels(), fontsize=16)\n",
    "    plt.setp(ax[i].get_yticklabels(), fontsize=16)\n",
    "\n",
    "ax[1].set_title(\"Our model\", fontsize=16)\n",
    "# ax[1].set_title(\"Sigma with prior\", fontsize=16)\n",
    "ax[0].set_title(\"NobBS\", fontsize=16)\n",
    "\n",
    "ax[0].set_ylabel(\"Absolute error\", fontsize=16)\n",
    "\n",
    "plt.suptitle(\"\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(f\"comparison_results.pdf\", dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "739aa322-6cf3-4fa5-813d-9c61fb6947ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAENCAYAAAASQ8TFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAydElEQVR4nO3de3QUZZ7/8U8n6VxIuASQAEaQIIRAYLwwwAgqCjIrXqLLxQEnihthAdFMECWwgyyCLiqIyGWQMFyMilyCBERMwPUIwhoW1JFLQATZDBjBQAIm5J7+/eGhf7YJdAerU5XO+3UO58hTT1d9O88hXz9d1VU2h8PhEAAAAADgN/MzuwAAAAAA8BUELAAAAAAwCAELAAAAAAxCwAIAAAAAgxCwAAAAAMAgBCwAAAAAMAgBC/BRo0ePVt++fRUdHW12KQAAuLVz507FxcUpOjpaCxYsqNVr58+fr7vvvlvR0dHKysryUoWAZwLMLgBoyI4fP66kpCTl5ubq4sWLevfdd9WjRw+XOc8884y+/vpr5eXlqV27dpo8ebJuvfVWt/tOSUnRggULtHDhwivOy87OVnJysnJyclReXq6OHTtKkn766SdJ0sCBA/X0008rLCzM+ZqqqiqlpaVp7dq1Ki8vV1VVlSorK9WpUyf169dPQ4cOre2PAgAg6fDhw0pJSdHhw4cVEBCgiooKde7cWaNHj1bXrl3NLs+rbrvtNt12221X9cFgYmKi+vTpo0cffdQLlQG1wxkswERRUVFKT0/XXXfdpfLyck2cOFGFhYUuc+bOnatZs2YpNjZW6enpHoWr2oiJiVF6erpiY2PVqlUrpaenKz09Xf/93/+t559/XqmpqXrxxRddXvPaa69p/vz5mj17tjZu3KhNmzZp1apVOnfunN544w1D6wOAhuKDDz5QfHy8BgwYoM2bNys9PV2bN2/WgAEDFB8fr/fff9/sEgF4gIAFWMSAAQP0z3/+U9OnTze7FKf+/furU6dO+uSTT1zG161bp8GDBzvPdklSy5YtlZycXNclAoBPOHz4sJKTk5WcnKzBgwfLz+/n/0Xz8/PTfffdp+TkZE2bNk0HDx40uVIA7hCwAIsYOHCg4uPj9cEHH2jDhg1u53/88ccaMmSIBg0apDvvvFNTpkzR2bNna5z79ddfKz4+XnfffbcGDhyojRs3elxXZWWls9FfUlFRodzc3GpzY2JitHbtWo/3DQD42eLFixUSEqIHHnigxu1xcXEKDg7WkiVLJEkrV66s9p2j3bt3a/DgwYqOjnbpI3FxcerVq5fuuusu7dy5U/Hx8c5L8S5cuFDtWKdPn1ZcXJxuuukmxcfHa8OGDRoxYoRuvfVW/fWvf1V5ebm2b9+u+Ph49evXT88//7zKyspc9lFQUKDnn39ed911l/74xz8qLi5OH374YbVjbdiwQQMHDtTdd9+t+Ph4ff311zW+/++//16JiYm68847nf3ycnMBsxGwAAt57rnn1K1bN82cOVPffffdZed9+OGHeuqppzRu3DhlZmbqo48+0pkzZxQfH6+LFy9Wm79q1SqlpKRo27ZtGjdunCZPnqxdu3ZdsZaqqiqtW7dOx44d0xNPPOGyrXfv3srMzNTUqVN1+PBh57jNZlPr1q1r+a4BoGGrqqrSzp071a1bN9nt9hrnBAYGqmvXrvrss89UVVWlUaNGadasWS5zbr31Vi1durTaay9dil5QUKAdO3borbfeUmZmplq2bFnjsSIiIpyXjh87dkwVFRVavXq1UlNTlZaWppkzZ+r8+fNKTU3VihUrtG7dOqWnpztfX1ZWpscff1zHjx/X5s2blZGRoWeffVbPPfec1q9f75y3Y8cOTZkyRePGjdO2bduUkpKilStXVqsnPz9fI0eOVHl5uTIzM7V9+3b1799fjz76qI4dO+bJjxioUwQswEICAwM1f/58+fn5aeLEidU+EZQkh8OhV155Rb1799bAgQMlSUFBQXr22Wd17NgxrVmzptpr/vznPys4OFiSNGTIEN1www013vzizJkziouLU1xcnP7whz9o1qxZGjdunB577DGXedOnT9fvf/97paWlKS4uTgMGDNBLL73kErYAAJ7Jz8/XxYsX1aJFiyvOa9mypS5evKj8/PyrOk5RUZHGjh0rm82mkJAQvffeey43MKpJVVWVhgwZIknq2LGjbrjhBmVmZurBBx+UJHXq1EkdO3bU559/7nxNenq6Dh06pKSkJIWGhkqS+vXrp/79+2vOnDnO3rZ48WJFRUU59x8cHKxHHnmkWg0rV65Ubm6uJk+e7Aygo0aNUpMmTZSSknJVPwvAmwhYgMVcd911mjVrlg4dOqRXX3212vbjx48rNzdX3bt3dxnv0qWLgoKCajwz1alTJ5e/d+vWTfv371dVVZXL+C9vcpGVlaUNGzbo448/1siRI1VSUuKcFxERobffflvr1q3TY489Jn9/f61atUpxcXHVPlEFABjLZrNd1euaNWvmEuKuu+66apeA/1pkZKT8/f2df2/atGm1sWbNmikvL8/59927d0tStT7Vo0cP5efnKzs7W5WVldq/f3+1OyN27ty5Wg3/8z//oxYtWqh9+/bOMX9/f0VFRWnPnj1XrB8wA7dpByzonnvuUVZWllJTU9W3b1+FhIQ4t1365LJJkybVXtekSZMaP9n89SeUTZo0UXl5uc6dO3fZS0Sknz+tnDhxov793/9d69atU3x8vMv2Hj16qEePHpo6dar27dunGTNmKDU1VXfeeaf69u1bq/cMAA1Vs2bN1KhRo8t+j/aSvLw8hYaGqlmzZld1nEtnk2rjl/1HkvPs16/HfvmBXX5+vkJCQhQYGOgyr2nTps7t+fn5qqioqNbLGjduXK2G/Px8FRYWKi4uzmX8p59+ksPhqPV7AryNM1iARU2ZMkXR0dFKTk7Wjz/+6BwPDw+XJJ0/f77aay5cuODc/ku/vvX7+fPnZbfb1bx5c7d1XHfddZKkb7/91jm2adMmVVZWusy75ZZb9J//+Z+SpEOHDrndLwDgZ/7+/urXr58OHjyo8vLyGueUlZXp0KFDuv32251nnS6dRfplyCguLvZ+wW6Eh4eruLi42mXul/pWeHi4wsPDZbfbq91ko6abboSHh+uaa65xXmHxy8eJ/Pout4AVELAAiwoKCtLrr7+u0tJS/dd//ZdzPCoqSm3atNH+/ftd5h85ckSlpaU1njk6evSoy98PHjyo7t27u700RJJ++OEHSXI50/X666/rxIkT1eZe2l9NIQ8AcHnjx49XcXGxNm3aVOP2jRs3qqSkRGPHjnWOXfqQ7Jeh5Eo3SKorl57X+Ou7/O3fv1/h4eGKiYmRv7+/YmNjq912/tf96tL+Tp8+rYKCApfxnTt36vXXXze0dsAIBCzAwjp06KAXXnjB5dp2m82m5557TllZWdq+fbskqbS0VHPmzFFUVJQefvjhavtJSUlxfocqLS1Nx44d04QJE9we/9y5c1q8eLEaNWrk/ELzJS+99JJLXWfOnNHcuXPVqlUrDRo06GreLgA0WDExMZo9e7Zmz56tDz/80HnJXVVVlbZs2aJXX31VL7/8srp06eJ8Tbt27dS6dWtnLygpKdHmzZtNqf+X4uLi1LVrV82bN09FRUWSfv5e1ieffKJnnnnGeenghAkT9N133yktLU3Sz/XXdNOKUaNGqVWrVnrppZecZ8Vyc3P14osvuvw8AKuwObh4FTDN8ePHlZSUpNzcXDVq1Eht27bVu+++W23eX//6V/3f//2fUlNTnWPbt2/X4sWLVVhYqLKyMvXp00fPPvus8wvMo0eP1qFDh5SXl6cVK1bojTfecAaiCRMmOANTdna2kpOTlZOTo/LycufDg8vLy1VSUqLu3btr3LhxLk1s+/bt2r59uw4ePCh/f38VFxeroqJCvXr10pNPPqnIyEhv/cgAwKdlZ2dr6dKlOnLkiOx2u8rLy9WlSxeNGTOmxjCRlZWlWbNmqbS0VFFRUXrsscc0atQotWnTRrfffrteeOEFxcfH68iRI7p48aI6duyo4cOH13i3vktKSkr08MMPKycnR5J0ww03aN26dRo2bJjzcvGaxtq1a6c1a9YoODhYBQUFmjt3rj777DMFBgYqODhYY8aM0b333utyrI0bNzrvatuyZUslJiZq1KhRatmypbp27eoMXD/88INeffVV7d27V+Hh4QoKCtJjjz2mwYMHS5Lmz5+vDz74QDk5OWrXrp3i4uI8+iAR8AYCFgAAAAAYhEsEAQAAAMAgBCwAAAAAMAgBCwAAAAAMQsACAAAAAIMEmHnwjIwMvfXWWyouLtb58+fVtGlTPfroo867m5WXl2vRokX66KOPFBAQoLCwME2aNEk9e/b0+BhfffWVgoKCvPQO6kZpaWm9fw/1HWtgDayDNfjCOpSWlurGG2/0aC59BEZhHayBdbAGX1iHy/US0wLWypUrlZ6err/97W9q3bq1ysvLlZycrM8//9wZsGbOnKmsrCytXr1azZs317p165SQkKD33ntPMTExHh0nKCjI47lWlZ2dXe/fQ33HGlgD62ANvrAO2dnZHs+lj8AorIM1sA7W4AvrcLleYsolgidPntScOXM0Y8YMtW7dWpJkt9v13HPPOZ/LcPz4ca1du1ajR492Pql82LBhioyM1Lx588woGwAAAACuyJSAlZ6eriZNmqhHjx4u4xEREerevbuknx9k6nA41Lt3b5c5ffr00a5du5xPBgcAAAAAqzDlEsEvv/xS1157rTIyMrRq1SqdO3dOzZo109ChQzV06FBJ0pEjR+Tn56e2bdu6vDYyMlIVFRU6duxYtYBWk9LS0lpdCmJFJSUl9f491HesgTWwDtbQ0NaBPgKjsA7WwDpYgy+vgykB64cfftCpU6e0fPlyLVy4UC1atFBGRoaeeeYZ/fjjjxo3bpzy8/MVHBwsf39/l9eGhYVJkvLz8z06FtfOwwisgTWwDtbgC+vAd7BgBtbBGlgHa/CFdbDUd7BKS0t18eJFPfvss7rmmmvk5+ene+65RwMGDNCbb76p4uJiM8oCAAAAgN/ElIAVGhoqSdVSa0xMjIqLi/Xtt98qPDxcJSUlqqysdJlTWFgoSQoPD6+bYgEAAADAQ6YErKioKEmSw+FwGb90OaDD4VB0dLSqqqqUm5vrMufkyZMKCAhQx44d66ZYAAAAAPCQKQHrzjvvlCQdPnzYZfybb75RcHCwbrjhBg0cOFA2m0179uxxmZOVlaW+ffs6z4LVV7GxsbLZbB796dq1q8dzY2NjzX5rAAAAQINlSsAaPHiwunfvrvnz5ztvt753715lZGRo7NixatSokaKiojR8+HC9+eabOnfunCQpLS1NOTk5SkpKMqNsQx04cEAOh8OjP+0nf+Dx3AMHDpj91gAAAIAGy5S7CPr7+2vZsmWaM2eO7rvvPgUGBiowMFDPP/+8hg8f7pw3bdo0LVy4UCNGjJDdbldoaKiWL19e7+84AgAAAMA3mRKwJKlZs2aaNWvWFefY7XYlJSX5xBkrAAAAAL7PlEsEAQAAAMAXEbAAAAAAwCAELAAAAAAwCAELAAAAAAxCwAIAAAAAgxCwAAAAAMAgBCwAAAAAMAgBCwAAAAAMQsACAAAAAIMQsAAAAADAIAQsAAAAADAIAQsAAAAADELAAgAAAACDELAAAAAAwCAELAAAAAAwCAELAAAAAAwSYNaBT548qfvvv1/t2rWrti01NVVNmjSRJJWXl2vRokX66KOPFBAQoLCwME2aNEk9e/as65IBAAAA4IpMC1iSFBsbq9TU1CvOmTlzprKysrR69Wo1b95c69atU0JCgt577z3FxMTUUaUAAAAA4J6lLxE8fvy41q5dq9GjR6t58+aSpGHDhikyMlLz5s0zuToAAAAAcGXpgLV9+3Y5HA717t3bZbxPnz7atWuXioqKTKoMAAAAAKoz9RLBvLw8TZo0SSdOnNBPP/2k2NhYjRkzRtHR0ZKkI0eOyM/PT23btnV5XWRkpCoqKnTs2DH16NHDjNIBAAAAoBrTApa/v78CAgI0atQoxcbGqrCwUDNnztTw4cOVmpqqHj16KD8/X8HBwfL393d5bVhYmCQpPz/f7XFKS0uVnZ3tlfdQl3zhPdRnJSUlrIEFsA7W0NDWwRf6SENbM6tiHayBdbAGX14H0wJWmzZttHnzZuffw8LCNGPGDH366aeaN2+eVqxYYchxgoKCfOBmGMd94D3Ub9nZ2ayBBbAO1uAL61Cbpu4LfcQX1swXsA7WwDpYgy+sw+V6iaW+gxUcHKzOnTvrq6++kiSFh4erpKRElZWVLvMKCwud2wEAAADAKkwLWD/99JPKysqqjfv7+6uqqkqSFB0draqqKuXm5rrMOXnypAICAtSxY8c6qRUAAAAAPGFawHrxxReVmZnpMlZWVqZvvvlGXbt2lSQNHDhQNptNe/bscZmXlZWlvn37KjQ0tM7qBQAAAAB3TL1E8O9//7vOnDkjSaqsrNQrr7yic+fOacKECZKkqKgoDR8+XG+++abOnTsnSUpLS1NOTo6SkpJMqxsAAAAAamLaTS4ef/xxrVmzRk888YSkn+8I2LFjR61YsUJ9+vRxzps2bZoWLlyoESNGyG63KzQ0VMuXL6/3X4oDAAAA4HtMC1jR0dF6/vnn3c6z2+1KSkrijBUAAAAAy7PUXQQBAAAAoD4jYAEAAACAQQhYAAAAAGAQAhYAAAAAGISABQAAAAAGIWABAAAAgEEIWAAAAABgEAIWAAAAABiEgAUAAAAABiFgAQAAAIBBCFgAAAAAYBACFgAAAAAYhIAFAAAAAAYhYAEAAACAQQhYAAAAAGAQAhYAAAAAGMQyAWvkyJGKjo7WyZMnzS4FAAAAAK5KgNkFSFJGRob27dtX47aioiLNnTtXu3btkr+/v1q3bq0pU6aoU6dOdVwlAAAAAFyZ6WewysrKNHfuXN1xxx01bk9MTFR2drY2btyoDz/8UD169FB8fLxOnz5dx5UCAAAAwJWZHrDeffddde/eXd27d6+2bdeuXdq5c6eefvpphYSESJLGjx+vqqoqLVmypK5LBQAAAIArMjVgFRQUaNmyZZo4cWKN2zMyMmS323XLLbc4xwIDA3XzzTcrMzOzrsoEAAAAAI+YGrAWLVqkBx54QNdee22N27/55hu1atVKgYGBLuORkZHKy8vT2bNn66JMAAAAAPCIaTe5OHHihLZu3aqtW7dedk5+fr5CQ0OrjV8aKygoUIsWLa54nNLSUmVnZ/+2Yi3AF95DfVZSUsIaWADrYA0NbR18oY80tDWzKtbBGlgHa/DldTAtYM2ZM0djxoxR48aNvXqcoKAgxcTEePUY3nfcB95D/Zadnc0aWADrYA2+sA61aeq+0Ed8Yc18AetgDayDNfjCOlyul5hyieDevXt19OhRjRgx4orzwsPDVVRUVG380lizZs28UR4AAAAAXBVTzmDt2rVLlZWVGjp0qHMsLy9PkjRmzBjZ7XZNnDhRnTt31oEDB1RWVubyPayTJ0+qZcuWbi8PBAAAAIC6ZErASkxMVGJiosvYggULtHDhQi1dulSRkZGSJH9/f61Zs0ZffvmlevfuLenn52Z98cUXuvfee+u8bgAAAAC4EtOfg3Ul/fr1U79+/TR//nwVFxdLkpYsWSI/Pz+NHTvW5OoAAAAAwJVpN7m45NNPP9Vrr71W7RLB9PR0SdIbb7yhOXPm6MEHH5Sfn59at26t1NRURUREmFk2AAAAAFRjesC64447dMcdd1x2e2hoqKZPn16HFQEAAADA1bH0JYIAAAAAUJ8QsAAAAADAIAQsAAAAADAIAQsAAAAADELAAgAAAACDuA1YOTk5Wr9+vdavX+98FpUkHTlyRCNHjtRNN92ku+66S8uXL/dqoQAAAABgdW4D1po1azR9+nR98sknqqiokCQVFxdr7Nix+vLLL3XzzTcrKipKr732mjZs2OD1ggEAAADAqtw+B2vfvn169dVXNXjwYOfY1q1blZubq8TERI0bN06StHnzZr3zzjv613/9V+9VCwAAAAAW5vYMVlFRkUu4kqQPP/xQwcHBio+Pd47df//9OnPmjPEVAgAAAEA94TZgBQS4nuS6ePGi9uzZoz59+igsLMxlW+PGjY2tDgAAAADqEbcBq6qqSpWVlc6/Z2Zmqry8XAMGDPBqYQAAAABQ37gNWDExMXr99ddVVFSko0ePasGCBQoJCal22WBmZqaaNGnitUIBAAAAwOrc3uRi3LhxGj58uJYtW+YcmzFjhkJDQyVJX3/9tVauXKnMzExNmDDBe5UCAAAAgMW5DVjt27fXBx98oIyMDFVUVKhXr17q2rWrc3vjxo1122236bbbbtOtt97q1WIBAAAAwMrcBixJuuaaa/TnP/+5xm0dOnRQhw4dDC0KAABcWWxsrA4ePGj4frt166YDBw4Yvl8AaCg8CljFxcXKzMxUXl6eOnTooP79+8vPz+3XtwAAgJfUJgRdn7xFJ2bf68VqAACXuA1YP/74o/70pz/p+++/l8PhkM1mU7du3bRy5cpqt2kHAAAAgIbMbcB67bXXZLfbNWPGDLVp00YnTpzQ8uXLtWTJEk2aNOmqDpqTk6PVq1crKytL0s8PM27RooXGjBmj/v37O+eVl5dr0aJF+uijjxQQEKCwsDBNmjRJPXv2vKrjAgAAAIA3uQ1Y+/btU2pqqiIiIiRJt912m3r16qX/+I//uOqD7tixQ1u2bFFqaqrat2+vqqoqzZkzR+PGjdOqVavUq1cvSdLMmTOVlZWl1atXq3nz5lq3bp0SEhL03nvvKSYm5qqPDwAAAADe4PaLVHa73RmuLomOjlZxcXG1ud9//71HB42IiNBTTz2l9u3b/1yEn5/Gjh2rqqoqffzxx5Kk48ePa+3atRo9erSaN28uSRo2bJgiIyM1b948j44DAAAAAHXJbcAKDAz0ePzJJ5/06KB33323hg0b5jJWWFgoSc4wtX37djkcDvXu3dtlXp8+fbRr1y4VFRV5dCwAAAAAqCtuLxEsKCjQxo0bq42fP3++2vj58+evqojTp0/rhRdeULdu3fTII49Iko4cOSI/Pz+1bdvWZW5kZKQqKip07Ngx9ejR46qOBwAAAADe4DZg5ebmasqUKXI4HNW2JScnu/zdZrPV6uA5OTlKSEhQTk6Obr/9di1atMh5Z8L8/HwFBwfL39/f5TW/3O6J0tJSZWdn16ouK/KF91CflZSUsAYWwDpYQ0NbB/oIjNLQ/u1YFetgDb68Dm4DVrt27TRr1iy3O3I4HJo2bVqtDt6uXTtt27ZNhYWFmj9/vuLi4rR48WJD7xIYFBTkAzfEOO4D76F+y87OZg0sgHWwBl9Yh9o0dfoIjOIL/3Z8AetgDb6wDpfrJW4D1h/+8AfnXf08mXs1wsLCNHXqVH322WeaMWOGNm/erPDwcJWUlKiystLlLNal72qFh4df1bEAAAAAwFvc3uRixowZHu/scjfE+LWSkpJqlxzabDZ17txZR48eVVlZmaKjo1VVVaXc3FyXeSdPnlRAQIA6duzocV0AAAAAUBfcBixPlZaWaseOHR7NHT16tL766qtq46dOnVJYWJjsdrsGDhwom82mPXv2uMzJyspS3759FRoaakTZAAAAAGAYjwLWsWPHtHDhQr344ovKzMx02VZYWKg333xTd911l3Jycjw+8IIFC5w3qnA4HHrrrbe0f/9+xcfHy2azKSoqSsOHD9ebb76pc+fOSZLS0tKUk5OjpKQkj48DAAAAAHXF7Xew/vd//1cJCQkqKyuTJL399tt6+umn9cQTT2jRokV6++23dfHiRd1111169NFHPTroX/7yF61fv17x8fHy9/dXaWmpmjVrpldffVX333+/c960adO0cOFCjRgxQna7XaGhoVq+fHm9/0IcAAAAAN/kNmDNnz9fgwcP1mOPPaby8nItWbJEb7/9tk6dOqWPPvpIw4YN0yOPPKLIyEiPD3rLLbfolltucTvPbrcrKSmJM1YAAAAA6gW3Aev777/XypUrFRDw89TZs2erd+/eKi8v18cff6ymTZt6vUgAAAAAqA/cBqzQ0FBnuJKkJk2aKCIiQrNmzZLdbvdqcQAAAABQn7i9ycUvw9Ul4eHhNYarf/u3fzOmKgAAAACoh9yewSoqKtLevXtdnlt18eLFamPSz7dZBwAAAICGym3AysnJUXx8vMuYw+GoccxmsxlbHQAAAADUI24DVps2bfT000+73ZHD4dDChQsNKQoAAAAA6iO3AatLly566KGHPNrZtm3bfnNBAAAAAFBfub3Jxd/+9rcrbi8qKlJVVZVHcwEAAADAl7kNWN9995169eqlXr16adq0adW2f/HFF/rjH/+ovXv3eqVAAAAAAKgv3Aasbdu2qaKiQpMnT9aUKVOqbe/Ro4cGDBigsWPH6tixY14pEgAAAADqA7cBa8eOHZo9e7aGDBmiRo0aVdvetGlTJScnKzExUSkpKV4pEgAAAADqA7cBKzc3V4MGDXK7oxEjRig7O9uQogAAAACgPnIbsMLCwjzaUUBAgPz83O4OAAAAAHyW20TkcDg82pHD4VBFRcVvLggAAAAA6iu3ASs6OloZGRlud7Rt2zZ16tTJkKIAAAAAoD5y+6DhhIQExcfHq6ioSHFxcfL393fZXllZqU2bNmn27Nl66623vFYoAAAAAFid24DVpUsXTZ06VdOmTdPLL7+srl27qkWLFpKks2fP6tChQyopKdHs2bMVHR3t9YIBAAAAwKrcBixJeuihh9SpUyctWbJEu3fv1sWLFyVJISEhuu222/Tkk0/WKlxlZ2frnXfe0b59+xQQEKDKykrdeuutGj9+vJo3b+6cV1RUpLlz52rXrl3y9/dX69atNWXKFC5FBAAAAGBJHgUsSYqNjdXChQtVVVWlgoICSVKzZs2u6s6BSUlJ6tSpk9LS0tSoUSOdPn1ajz32mHbu3Kn09HQFBwdLkhITE1VUVKSNGzcqJCREr7/+uuLj45Wenq6IiIhaHxcAAAAAvKnW6cjPz0/NmzdX8+bNf9Nt2SdNmuR8cHFERIQSEhJ04sQJffrpp5KkXbt2aefOnXr66acVEhIiSRo/fryqqqq0ZMmSqz4uAAAAAHiLKQ+u2rRpk9q3b+8y1qpVK0nShQsXJEkZGRmy2+265ZZbnHMCAwN18803KzMzs+6KBQAAAAAPeXyJoJECAwOrjZ04cUI2m009e/aUJH3zzTdq1apVtbmRkZH65JNPdPbsWefNNq6ktLRU2dnZxhRuIl94D/VZSUkJa2ABrIM1NLR1oI/AKA3t345VsQ7W4MvrYErA+rXKykqtX79eQ4cOVYcOHSRJ+fn5Cg0NrTb30lhBQYFHASsoKEgxMTHGFlznjvvAe6jfsrOzWQMLYB2swRfWoTZNnT4Co/jCvx1fwDpYgy+sw+V6iSmXCP7a4sWLZbfbNXXqVLNLAQAA8FhsbKxsNptHf7p27erx3NjYWLPfGoCrZHrASktL09atW5WSkuK86YUkhYeHq6ioqNr8S2PNmjWrqxIBAABqdODAATkcDo/+tJ/8gcdzDxw4YPZbA3CVTA1YGzdu1IoVK7Rq1apql/t17txZZ86cUVlZmcv4yZMn1bJlS48uDwQAAACAumRawEpPT9eyZcu0YsUKXXPNNZKkTz75RGvWrJEkDRo0SOXl5fryyy+drykrK9MXX3yhQYMGmVIzfAuXdQAAAMBoptzkYtOmTZo2bZoSExO1e/du5/i+ffucYatfv37q16+f5s+fr7///e8KCQnRkiVL5Ofnp7Fjx5pRNnxMbS6/uD55i07MvteL1QAAAMAXmBKwZs2apdLSUr3yyivVtk2YMMH532+88YbmzJmjBx98UH5+fmrdurVSU1MVERFRl+V65HczMnW+uNxr+78+eYvh+2waYtc/pnM2EAAAADCKKQFrz549Hs0LDQ3V9OnTvVyNMc4Xl3vtDIe3bmPpjdAGAAAANGSm30UQAAAAAHwFAQsAAAAADELAAgAAAACDELAAAAAAwCAELAAAAAAwCAELAAAAAAxCwAIAAAAAgxCwAAAAAMAgBCwAAAAAMAgBCwAAAAAMQsACAAAAAIMQsAAAAADAIAQsAAAAADAIAQsAAAAADELAAgAAAACDELAAAAAAwCCmBqwzZ84oISFB0dHRZpYBAAAAAIYIMOvAmZmZmj17tgICrlxCeXm5Fi1apI8++kgBAQEKCwvTpEmT1LNnzzqqFAAAAAA8Y9oZrJSUFC1fvlw333zzFefNnDlTW7du1bvvvqsPPvhAQ4YMUUJCgrKzs+uoUgAAAADwjGlnsFavXu327NXx48e1du1azZo1S82bN5ckDRs2TCtXrtS8efO0dOnSuigVAIA68bsZmTpfXO6VfV+fvMXwfTYNsesf0wcZvl8AqM9MC1juwpUkbd++XQ6HQ71793YZ79Onj9577z0VFRUpNDTUWyUCAFCnzheX68Tsew3fb3Z2tmJiYgzfrzdCGwDUd6YFLE8cOXJEfn5+atu2rct4ZGSkKioqdOzYMfXo0eOK+ygtLa2zywm9dZySkhKv7ZtLLT3Hz8p83vy3AM81tHWoyz4ieed3DX3EOvh5ma+h/Q6zKl9eB0sHrPz8fAUHB8vf399lPCwszLndnaCgIK98alfdca8dx1ufPHqzZt/Dz8oKvPdvAbXhC+tQm6Zed31E8tbvGvqIVfDzsgJf+B3mC3xhHS7XS3gOFgAAAAAYxNIBKzw8XCUlJaqsrHQZLywsdG4HAAAAAKuwdMCKjo5WVVWVcnNzXcZPnjypgIAAdezY0aTKAAAAAKA6SwesgQMHymazac+ePS7jWVlZ6tu3L3cQBAAAAGAplg5YUVFRGj58uN58802dO3dOkpSWlqacnBwlJSWZXB0AAAAAuDLtLoIvv/yydu/e7bz8Ly4uTpK0bt06BQYGOudNmzZNCxcu1IgRI2S32xUaGqrly5fX+7uOAPhZbGysDh48aPh+u3XrpgMHDhi+XwCA9dBLYCWmBazJkyd7NM9utyspKYkzVoCPqk3juj55i1cewgoAqN/oJbASSz8HC6it383I1Pnicq/s+/rkLYbvs2mIXf+YPsjw/QIAAMAcBCz4lPPF5V75VMpbD8PzRmgDAACAeQhYBmnUYZ66r0r23gH2uJ9SW406REjiFDkAAL/E1RAAfgsClkEufpfktet5OXsCAA2DVz+s44M6j3E1BIDfgoAFAIBFeOvDOv7HHgDqjqWfgwUAAAAA9QkBCwAAAMBvFhsbK5vN5tGfrl27ejw3NjbW7LdWKwQsAAAAAL/ZgQMH5HA4PPrTfvIHHs+tbw97JmABAAAAgEEIWAAAAABgEAIWAIDr5gEAMAgBCwDAdfMAABiE52AZyLvPAzlu+B6bhtgN3ycA4LfxXi+hjwBAXSBgGcQbD4a85PrkLV7dvy9p1GGeuq9K9s7O9xi/y0YdIiSxtgB+5q3f9fSR2qGXWMPvZmTqfHG5V/btjQ8ymobY9Y/pgwzfL+ofAhZ8ysXvkrzyPxHZ2dmKiYkxfL/ePetpHpoigPqMXmIN54vLWQfUSwQsAIajKQIAgIaKgGWS2NhYHTx40OP5tpc9m9etWze+VA4ADQB9xLv4LhyAq2X5gHX27Fm99NJLzl/20dHRmjp1qlq3bm1yZb9NbZqXtz6191U0RQANAX3Ee/guHPD/cdl/7Vk6YJWVlenxxx/X9ddfry1btshms2nq1Kl69NFH9f777ys0NNTsEmExNEVr4Avi1kBTBFCf0Uusgcv+a8/SAWvjxo06cuSIFi1apICAn0udNGmSbr/9dq1evVpPPPGEyRWiPuPyGu/hC+LWQFMEvI9e4j0Xv0syu4Ra8dWrUgi6tWfpgJWRkaG2bdvquuuuc45dc801uuGGG5SZmUnAwm/C5TXe5en/LH//9/Eqz8sx/Pj2lu3UNmGxR3NpilfBC02xccfWMrspArVFL/EerkqxBj40rT2bw+FwmF3E5fTr108dOnRQamqqy/jYsWP1+eef66uvvnK7j6+++kpBQUFeqrBulJSUKDg42OwyGjTWwBpYB+954IEH9O233xq+3xtuuEGbNm0yfL9G8LSx00dgFNbBexri77C6cs8qz7+DboUPTcMC/bRuxPWG13A5NfUSS5/BKigoqPF7VmFhYSouLvboF1VQUFC9/7SIT7zMxxpYA+vgPUePHvV4ri+sQ3Z2tsdz6SMwCuvgPQ3td1hdOjG7Fj+r2f/n8VRfWIfL9RK/Oq4DAAAAAHyWpQNWs2bNVFRUVG28sLBQISEhnGYHAAAAYCmWDljR0dE6depUtfFTp06pc+fOJlQEAAAAAJdn6YA1aNAgnTp1SidPnnSO5eXl6dixYxo0iGelAAAAALAWSweshx56SJ07d9acOXNUUVGhqqoqzZ07V5GRkRoxYoTZ5QEAAACAC0sHrMDAQK1YsUL+/v669957dc8996iwsFCrVq2q8e6CAAAAAGAmS9+mXZJatmypuXPnml0GAAAAALhl6QcNG8EXHhAJADBWaWmpbrzxRo/m0kcAADW5XC/x+YAFAAAAAHXF0t/BAgAAAID6hIAFAAAAAAYhYAEAAACAQQhYAAAAAGAQAhYAAAAAGISABQAAAAAGIWBZ1JkzZ5SQkKDo6GizSwGAakaOHKno6GidPHnS7FJwBfQSAFbly30kwOwCUF1mZqZmz56tgACWxyzZ2dl65513tG/fPgUEBKiyslK33nqrxo8fr+bNm5tdXoORk5Oj1atXKysrS5JUVFSkFi1aaMyYMerfv7+5xTVgGRkZ2rdvn9llwA16ifnoJdZAL7EeX+8jnMGyoJSUFC1fvlw333yz2aU0WElJSTp//rzS0tK0efNmrVixQp999plGjBihkpISs8trMHbs2KEtW7Zo3rx52rBhg7Zu3aobb7xR48aN0549e8wur0EqKyvT3Llzdccdd5hdCtygl5iPXmIN9BJraQh9hIBlQatXr9b1119vdhkN3qRJk9SoUSNJUkREhBISEnTixAl9+umnJlfWcEREROipp55S+/btJUl+fn4aO3asqqqq9PHHH5tcXcP07rvvqnv37urevbvZpcANeok10EvMRy+xlobQR7huwIK4nMN8mzZtUmBgoMtYq1atJEkXLlwwo6QG6e677642VlhYKElcXmOCgoICLVu2TGvWrNGGDRvMLgdu0EvMRy+xBnqJdTSUPsIZLKAGv26IknTixAnZbDb17NnThIogSadPn9YLL7ygbt266ZFHHjG7nAZn0aJFeuCBB3TttdeaXQpQL9BLrIleYp6G0kcIWIAHKisrtX79eg0dOlQdOnQwu5wGJycnR3fffbduv/12VVZWatGiRQoLCzO7rAblxIkT2rp1q8aNG2d2KUC9RS8xF73EXA2pjxCwAA8sXrxYdrtdU6dONbuUBqldu3batm2b9u3bp+uvv15xcXHau3ev2WU1KHPmzNGYMWPUuHFjs0sB6i16ibnoJeZqSH2EgAW4kZaWpq1btyolJcX5RWWYIywsTFOnTlWLFi00Y8YMs8tpMPbu3aujR49qxIgRZpcC1Fv0Euugl9S9htZH+AYscAUbN27UihUrtGrVKrVo0cLschqckpISBQUFyWazOcdsNps6d+6sjIwMlZWV1fgdBxhr165dqqys1NChQ51jeXl5kqQxY8bIbrdr4sSJPn3LXeC3oJeYi15ivobWRwhYwGWkp6dr2bJlWrFiha655hpJ0ieffKIzZ87o4YcfNrm6hmH06NGaOHGibrrpJpfxU6dOKSwsTHa73aTKGpbExEQlJia6jC1YsEALFy7U0qVLFRkZaVJlgPXRS8xHLzFfQ+sjBCygBps2bdK0adOUmJio3bt3O8f37dvnbJCoGwsWLNDcuXMVHh4uh8Oh1NRU7d+/X+PHj3f5NBIArIZeYh30EtQlm8PhcJhdBFy9/PLL2r17t3Jzc3X+/Hl16dJFkrRu3TpOYdeRXr166fz58zVumzBhgp566qk6rqhh2rdvn9avX6/9+/fL399fpaWlatasmUaOHKn777+fpmiCTz/9VK+99pry8vKUl5enjh07ym63Kz093ezS8Cv0EvPRS6yBXmItDaGPELAAAAAAwCDcRRAAAAAADELAAgAAAACDELAAAAAAwCAELAAAAAAwCAELAAAAAAxCwAIAAAAAgxCwAAAAAMAgBCwAAAAAMAgBC6hnzp49q7i4OPXq1UvR0dGKi4tTWlqa2WUBAOoJ+gjgXTaHw+EwuwgAtZecnKz3339fR44cMbsUAEA9RB8BvIMzWAAAAABgkACzCwDgPUuWLNH27dtVVVWl8vJyderUSZMnT1ZERITLvA0bNmjx4sWy2Wxq2bKlnn76aY0aNUotW7ZU165dlZKSYtI7AACYiT4C1B4BC/Bhy5Yt09tvv60uXbqosrJSL7/8ssaOHav169fL399fkrRjxw5NmTJFL730koYMGaKSkhL95S9/kST96U9/0lNPPWXiOwAAmIk+AtQelwgCPmzNmjXq0qWLJMnf318PP/ywDh06pAMHDjjnLF68WFFRURoyZIgkKTg4WKNHjzalXgCAtdBHgNrjDBbgwy5cuKDx48crJydH/v7+Ki8vlyT985//1O9+9ztVVlZq//79+pd/+ReX13Xu3NmMcgEAFkMfAWqPM1iAjzp8+LAeffRRtWvXThs3blR6erqWLl0qSSorK5Mk5efnq6KiQk2aNHF5bePGjeu8XgCAtdBHgKtDwAJ8UFlZmT788EOVlZVp7NixCgio+WR1eHi47Ha7Lly44DL+678DABoW+ghw9bhEEPBBCQkJKioqkiTZbDbneF5enss8f39/xcbG6uDBgy7jR48e9X6RAADLoo8AV48zWICPevbZZyVJK1askPTzp5E13SZ3woQJ+u6775SWliZJKikp0TvvvFN3hQIALIk+Alwdm8PhcJhdBADP/fjjj3r44YeVn5+vixcv6tprr61xzrJly3TixAktW7ZMfn5+at26tW6//Xa98soratOmje677z5NmjRJkrRx40YtXLhQNptNrVu31rPPPqthw4bpqaee0oQJE+r6LQIAvIg+AngXAQtANQUFBerdu7emT5+ukSNHml0OAKCeoY+gIeMSQaCB27t3r1avXu0yduna+U6dOplREgCgHqGPAK4IWEADV1BQoJSUFJ05c0aSVFhYqDfeeEM9e/bU73//e5OrAwBYHX0EcMVdBIEGLjo6WjfffLPi4+MVGBioCxcuqE+fPpo8ebLZpQEA6gH6COCK72ABAAAAgEG4RBAAAAAADELAAgAAAACDELAAAAAAwCAELAAAAAAwCAELAAAAAAzy/wDUN+ABcbez3wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(12, 4), sharey=True)\n",
    "\n",
    "pd.DataFrame(np.flip(np.array(our_crps_all))).boxplot(showfliers=False, whis=(5,95), ax=ax[1])\n",
    "pd.DataFrame(np.flip(np.array(nobbs_crps_all))).boxplot(showfliers=False, whis=(5,95), ax=ax[0])\n",
    "\n",
    "for i in range(2):\n",
    "\n",
    "    ax[i].set_xticklabels([1,2,3,4])\n",
    "    ax[i].set_xlabel(\"Lag\", fontsize=16)\n",
    "    plt.setp(ax[i].get_xticklabels(), fontsize=16)\n",
    "    plt.setp(ax[i].get_yticklabels(), fontsize=16)\n",
    "\n",
    "ax[1].set_title(\"Our model\", fontsize=16)\n",
    "# ax[1].set_title(\"Sigma with prior\", fontsize=16)\n",
    "ax[0].set_title(\"NobBS\", fontsize=16)\n",
    "\n",
    "ax[0].set_ylabel(\"CRPS\", fontsize=16)\n",
    "\n",
    "plt.suptitle(\"\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(f\"comparison_results_crps.pdf\", dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbf040ef-d51b-4b9c-a841-6b36c66f6f24",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}