docs/tutorials/wbic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# WBIC/BIC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook gives a tutorial on how to use Watanabe-Bayesian information criterion (WBIC) and Bayesian information criterion (BIC) for feature selection (Watanabe[2010], McElreath[2015], and Vehtari[2016]). The WBIC or BIC is an information criterion. Similar to other criteria (AIC, DIC), the WBIC/BIC endeavors to find the most parsimonious model, i.e., the model that balances fit with complexity. In other words a model (or set of features) that optimizes WBIC/BIC should neither over nor under fit the available data. \n",
    "\n",
    "In this tutorial a data set is simulated using the damped linear trend (DLT) model. This data set is then used to fit DLT models with varying number of features as well as a global local trend model (GLT), and a Error-Trend-Seasonal (ETS) model. The WBIC/BIC criteria is then show to find the true model. \n",
    "\n",
    "Note that we recommend the use of WBIC for full Bayesian and SVI estimators and BIC for MAP estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-28T21:49:16.127581Z",
     "start_time": "2022-04-28T21:49:13.501781Z"
    }
   },
   "outputs": [],
   "source": [
    "from datetime import timedelta\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import orbit\n",
    "from orbit.models import DLT,ETS, KTRLite, LGT\n",
    "from orbit.utils.simulation import make_trend, make_regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-28T21:49:16.132504Z",
     "start_time": "2022-04-28T21:49:16.129880Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1.3dev\n"
     ]
    }
   ],
   "source": [
    "print(orbit.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "## Data Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This block of code creates random data set (365 observations with 10 features) assuming a DLT model. Of the 10 features 5 are effective regressors; i.e., they are used in the true model to create the data.\n",
    "\n",
    "As an exercise left to the user once you have run the code once try changing the `NUM_OF_EFFECTIVE_REGRESSORS` (line 2), the `SERIES_LEN` (line 3), and the `SEED` (line 4) to see how it effects the results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-03-25T01:01:15.716689Z",
     "start_time": "2022-03-25T01:01:15.701163Z"
    }
   },
   "outputs": [],
   "source": [
    "NUM_OF_REGRESSORS = 10\n",
    "NUM_OF_EFFECTIVE_REGRESSORS = 4\n",
    "SERIES_LEN = 365\n",
    "SEED = 1\n",
    "# sample some coefficients\n",
    "COEFS = np.random.default_rng(SEED).uniform(-1, 1, NUM_OF_EFFECTIVE_REGRESSORS)\n",
    "trend = make_trend(SERIES_LEN, rw_loc=0.01, rw_scale=0.1)\n",
    "x, regression, coefs = make_regression(series_len=SERIES_LEN, coefs=COEFS)\n",
    "\n",
    "# combine trend and the regression\n",
    "y = trend + regression\n",
    "y = y - y.min()\n",
    "\n",
    "\n",
    "x_extra = np.random.normal(0, 1, (SERIES_LEN, NUM_OF_REGRESSORS - NUM_OF_EFFECTIVE_REGRESSORS))\n",
    "x = np.concatenate([x, x_extra], axis=-1)\n",
    "\n",
    "x_cols = [f\"x{x}\" for x in range(1, NUM_OF_REGRESSORS + 1)]\n",
    "response_col = \"y\"\n",
    "dt_col = \"date\"\n",
    "obs_matrix = np.concatenate([y.reshape(-1, 1), x], axis=1)\n",
    "# make a data frame for orbit inputs\n",
    "df = pd.DataFrame(obs_matrix, columns=[response_col] + x_cols)\n",
    "# make some dummy date stamp\n",
    "dt = pd.date_range(start='2016-01-04', periods=SERIES_LEN, freq=\"1W\")\n",
    "df['date'] = dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-03-25T01:01:15.730128Z",
     "start_time": "2022-03-25T01:01:15.720117Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(365, 12)\n",
      "          y        x1        x2        x3        x4        x5        x6  \\\n",
      "0  4.426242  0.172792  0.000000  0.165219 -0.000000 -0.171564  0.646560   \n",
      "1  5.580432  0.452678  0.223187 -0.000000  0.290559  1.760342  1.809456   \n",
      "2  5.031773  0.182286  0.147066  0.014211  0.273356 -1.536165 -1.029521   \n",
      "3  3.264027 -0.368227 -0.081455 -0.241060  0.299423  0.591493  0.505090   \n",
      "4  5.246511  0.019861 -0.146228 -0.390954 -0.128596 -1.447045 -1.678974   \n",
      "\n",
      "         x7        x8        x9       x10       date  \n",
      "0 -0.140274  2.024782 -0.334100 -1.179794 2016-01-10  \n",
      "1 -1.191844 -0.990802 -2.894241  1.757803 2016-01-17  \n",
      "2  0.572564 -0.750318  0.168999  0.524423 2016-01-24  \n",
      "3  1.348104  0.464408  1.263463 -1.705559 2016-01-31  \n",
      "4  0.292021  0.534810  1.463977 -0.137525 2016-02-07  \n"
     ]
    }
   ],
   "source": [
    "print(df.shape)\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WBIC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we use DLT model as an example. Different DLT models (the number of features used changes) are fitted and their WBIC values are calculated respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-03-25T01:05:53.061817Z",
     "start_time": "2022-03-25T01:01:15.732651Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n",
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 1 regressors: 1202.366\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 2 regressors: 1150.291\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 3 regressors: 1104.491\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 4 regressors: 1054.450\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 5 regressors: 1060.473\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 6 regressors: 1066.991\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 7 regressors: 1073.808\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 8 regressors: 1080.625\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 1000 and samples(per chain): 1000.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 9 regressors: 1084.251\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value with 10 regressors: 1088.111\n",
      "------------------------------------------------------------------\n",
      "CPU times: user 3.12 s, sys: 4.08 s, total: 7.2 s\n",
      "Wall time: 4min 37s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "WBIC_ls = []\n",
    "for k in range(1, NUM_OF_REGRESSORS + 1):\n",
    "    regressor_col = x_cols[:k]\n",
    "    dlt_mod = DLT(\n",
    "        response_col=response_col,\n",
    "        date_col=dt_col,\n",
    "        regressor_col=regressor_col,\n",
    "        seed=2022,\n",
    "        # fixing the smoothing parameters to learn regression coefficients more effectively\n",
    "        level_sm_input=0.01,\n",
    "        slope_sm_input=0.01,\n",
    "        num_warmup=4000,\n",
    "        num_sample=4000,\n",
    "        )\n",
    "    WBIC_temp = dlt_mod.fit_wbic(df=df) \n",
    "    print(\"WBIC value with {:d} regressors: {:.3f}\".format(k, WBIC_temp))\n",
    "    print('------------------------------------------------------------------')\n",
    "    WBIC_ls.append(WBIC_temp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also interesting to see if WBIC can distinguish between model types; not just do feature selection for a given type of model. To that end the next block fits an LGT and ETS model to the data; the WBIC values for both models are then calculated.  \n",
    "\n",
    "Note that WBIC is supported for both  the 'stan-mcmc' and 'pyro-svi' estimators. Currently only the LGT model has both. Thus WBIC is calculated for LGT for both estimators. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-03-25T01:06:50.978628Z",
     "start_time": "2022-03-25T01:05:53.065054Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 225 and samples(per chain): 25.\n",
      "WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
      "To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
      "INFO:orbit:Using SVI (Pyro) with steps: 301, samples: 100, learning rate: 0.1, learning_rate_total_decay: 1.0 and particles: 100.\n",
      "INFO:root:Guessed max_plate_nesting = 2\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value for LGT model (stan MCMC): 1144.606\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:orbit:step    0 loss = 303.3, scale = 0.11508\n",
      "INFO:orbit:step  100 loss = 116.73, scale = 0.50351\n",
      "INFO:orbit:step  200 loss = 116.58, scale = 0.5112\n",
      "INFO:orbit:step  300 loss = 116.84, scale = 0.50453\n",
      "INFO:orbit:Sampling (PyStan) with chains: 4, cores: 8, temperature: 5.900, warmups (per chain): 225 and samples(per chain): 25.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value for LGT model (pyro SVI): 1124.428\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated\n",
      "WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WBIC value for ETS model: 1203.031\n",
      "CPU times: user 50.9 s, sys: 8.06 s, total: 59 s\n",
      "Wall time: 57.9 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "lgt = LGT(response_col=response_col,\n",
    "          date_col=dt_col,\n",
    "          regressor_col=regressor_col,\n",
    "          seasonality=52,\n",
    "          estimator='stan-mcmc',\n",
    "          seed=8888)\n",
    "WBIC_lgt_mcmc = lgt.fit_wbic(df=df) \n",
    "print(\"WBIC value for LGT model (stan MCMC): {:.3f}\".format(WBIC_lgt_mcmc))\n",
    "\n",
    "lgt = LGT(response_col=response_col,\n",
    "          date_col=dt_col,\n",
    "          regressor_col=regressor_col,\n",
    "          seasonality=52,\n",
    "          estimator='pyro-svi',\n",
    "          seed=8888)\n",
    "WBIC_lgt_pyro = lgt.fit_wbic(df=df) \n",
    "print(\"WBIC value for LGT model (pyro SVI): {:.3f}\".format(WBIC_lgt_pyro))\n",
    "\n",
    "ets = ETS(\n",
    "        response_col=response_col,\n",
    "        date_col=dt_col,\n",
    "        seed=2020,\n",
    "        # fixing the smoothing parameters to learn regression coefficients more effectively\n",
    "        level_sm_input=0.01,\n",
    "    )\n",
    "\n",
    "WBIC_ets = ets.fit_wbic(df=df) \n",
    "print(\"WBIC value for ETS model: {:.3f}\".format(WBIC_ets))\n",
    "\n",
    "WBIC_ls.append(WBIC_lgt_mcmc)\n",
    "WBIC_ls.append(WBIC_lgt_pyro)\n",
    "WBIC_ls.append(WBIC_ets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot below shows the WBIC vs the number of features / model type (blue line). The true model is indicated by the vertical red line. The horizontal gray line shows the minimum (optimal) value. The minimum is at the true value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-03-25T01:06:51.187435Z",
     "start_time": "2022-03-25T01:06:50.982429Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyoAAAGbCAYAAAArwGQ+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAxOAAAMTgF/d4wjAAAwbklEQVR4nO3de7hcdX33/fd3Z+cAbgRJgwEjRZsCGoMhoIgPyqF6R0Hv0nArKmo9E6E3RWofxcdSrG1RW/HQ4k1QLKL0EYWIB5RUIShUQSUGY1QwYowB0lRQJAgkO/t7/zFrwmRnH2fPYc3M+3Vdc2XPOn5/syd71mfW+q1fZCaSJEmSVCZ97S5AkiRJkoYzqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqnf52F9AqM2fOzDlz5rS7DKm8hoZg69bHng8MQJ/fZUiSpOa4++67t2XmzNHm90xQmTNnDps2bWp3GVJ5bdkCH/vYY8/POAP226999UiSpK4WEf891ny/LpUkSZJUOgYVSZIkSaVjUJEkSZJUOj3TR0WSJElqlaGhITKz3WW0XUTQV+fNeQwqkiRJUoNs27aNjRs3sn379naXUhrTp0/nwAMPZMaMGZNaz6AiSZIkNcjGjRvZa6+9mD17NhHR7nLaLjO577772LhxI/Pnz5/UugYVSZIkqQGGhobYvn07s2fPpr/fw+yq2bNnc//99zM0NDSpy8DsTC9JkiQ1QLVPimdSdlV9PSbbZ8egIkmSJKl0DCqSJElSFzv//PN55JFH6lr3wx/+MJs3b95lW2effXaDKhubQUWSJEnqYu95z3tGDCqDg4Pjrjs8qLSSvXwkSZKkZhgchPvvb/5+9t0XRum8v2zZMgCe97znMW3aNA444ADmzp3L+vXr2bJlCz/96U+JCH7zm9+wzz77APAHf/AHfP/73+fyyy/nnnvu4dRTT2WPPfbgsssuA+Dee+/lpS99KT//+c+ZO3cuV111Ffvuu2/Dm2VQkSRJkprh/vvhYx9r/n7OOAP222/EWRdffDHLly/npptuYp999uF1r3sdt912GzfffDN77bXXmJs977zz+OQnP8mVV17JokWLALjmmmu49dZbue2225g9ezaveMUrWL58Oeeee26jW+WlX5IkSVIvednLXjZuSBnLi170ImbPng3A0Ucfzc9//vNGlbYLg4okSZLUQwYGBnZ5Pm3aNHbs2LHz+Xgd72fNmrXLuhPp61KPpl/6FREfBf4n8IfA4Zm5JiJmAZ8Fng48DGwB3pqZ64t1biyWf6DYzKcy80PFvP2Ay4E/Ah4FzsjMbzW7HZIkSdKk7Ltv5bKsVuxnDHvttRcPPPDAzj4ow82fP59bb72VE088kRUrVvDQQw/tnPf4xz+eBx54YMT1mq0VfVSuAj4A3Dxs+iXA1zIzI+IvgE8Ax9XMf1tmXjPC9t4H3JKZL4qIZwFfiIinZOb2xpcuSZIk1am/f9S+I630V3/1V7zwhS9kzz335IADDtht/oc+9CHOOuss3v3ud3PSSSftvKwL4KyzzuLNb34ze+65587O9K0Skx0hsu4dRWwATs7MNSPMOxK4KjMPKp7fCHx4pKASEVuB+Zm5uXj+XeBdmfmNsfY/b9683LRp09QaIXWzLVt27fA3Rsc8SZK0ux07dnDnnXdy8MEHM23atHaXUxqjvS4RcXdmzhttvbL0UflL4IvDpr0vItZGxJUR8VSAiJgNTK+GlMIG4MDhG4yIcyJiU/WxdevWZtUuSZIkTdnKdZtZ/N6vs3Jde8YtKZu2B5WIeBcwH6i9p9lrMvNQ4DDgJuArk91uZl6YmfOqj+GdhiRJkqSyWLluM2desZr7H9rGmVesNqzQ5qASEW8HlgIvzszfV6dn5q+KfzMz/xV4akTMzsz7gMGImFuzmYOAjS0sW5IkSWqYakgZHKp0yRgcSsMKbQwqEXEO8ErghZn525rp/RHxxJrnpwD/VYQUgM8Dy4p5zwKeBHyzVXVLkiRJjTI8pFQZVloQVCJieURsAuYBKyNifUTMAz4I7AOsiog1EXFrscpM4Nqif8rtwBlUbm9c9Q7guRHxM+Ay4NXe8UuSJEmdZrSQUtXrYaXptyfOzNNHmRWjLP8QcOQY2/sv4H80oDRJkiSpLcYLKVXVsHLRaYtZsmDumMt2m7Z3ppckSZJ6yURDSlWvnlkxqEiSJEktMtmQUlVvWHnPe97Dm970pp3Pb775ZiKCG2+8cee0ZcuW8Td/8zccdNBBHHLIISxatIinPe1pvOpVr9o5Sv2NN97IokWLdq6zdetWzj77bObPn88zn/lMDj/8cN7+9rezfXvjemQYVCRJkqQWqDekVNUTVo4//vhdQsmqVas46qijdpt2wgknAHDllVeyZs0a1q1bxwMPPDDiaPSZyUte8hIeeugh1q5dy+233853v/td5s+fz6OPPlpX20ZiUJEkSZJa4NwVa+sOKVWDQ8m5K9ZOePnnPOc53HPPPWzatAmonBk577zzdgaVe++9l40bN3L00Ufvst62bdv4/e9/zxOe8ITdtnnDDTewfv16LrroIvbYYw8Apk+fzrJly2jk2IUGFUmSJKkFLli6kP6+Ee8nNWH9fcEFSxdOePkZM2bw3Oc+l1WrVvHoo4/yi1/8ghNPPJFNmzbxyCOPsGrVKo4++mhmzZoFwKmnnsqiRYuYO3cufX19vPzlL99tm7fddhtHHHEEM2bMmFJbxmNQkSRJklpgyYK5XHTa4rrDSn9f1HX3r+rlX7feeivPfvazgcqZlu985zvceOONHH/88TuXrV769etf/5qDDjqId7zjHXXV2ggGFUmSJKlF6g0r9YYUqASVVatWsWrVKo477jgAjj322J3Tqv1Tdtlffz+nnHIK11133W7zjjjiCFavXs22bdsmXctkGFQkSZKkFppsWJlKSAF41rOexZYtW7jiiit2CSqf/exnuffee3eeZRnuhhtu4JBDDtlt+gknnMBTnvIUzjrrLB555BEABgcHueSSS9i6dWtdNY7EoCJJkiS12ETDylRDClQ6uh9zzDE8+OCDHHrooQAcfPDBPPjggxxzzDFMnz5957LVPirPeMYz+MlPfsJHPvKR3bYXEVx77bXMmDGDBQsW8IxnPIPDDjuMn/3sZzv7ujRCZE7tzgOdYt68eVm924GkEWzZAh/72GPPzzgD9tuvffVIktRhduzYwZ133snBBx/MtGnTJrTOWLcsbkRIKYPRXpeIuDsz5422nmdUJEmSpDYZ7cxKt4SUqTCoSJIkSW00PKwYUioMKpIkSVIDRFSCRj1dK6phZd/Hzei6kFJ9Paqvz0T1N6MYSZIkqdf09fUxffp07rvvPmbPnj3pA/MXHDqHF7yrcqvgHTt2NKPElstM7rvvPqZPn05f3+TOkRhUJEmSpAY58MAD2bhxI/fff3+7SymN6dOnc+CBB056PYOKJEmS1CAzZsxg/vz5DA0N1XUJWLeJiEmfSakyqEiSJEkNVu/BuR7jKyhJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiSpdAwqkiRJkkrHoCJJkiT1kJXrNrP4vV9n5brN7S5lTAaVFuqUN4UkSZK608p1mznzitXc/9A2zrxidamPSw0qLdJJbwpJkiR1n+rx6OBQAjA4lKU+LjWotECnvSkkSZLUXYYfj1aV+bjUoNJknfimkCRJUvcY7Xi0qqzHpQaVJurUN4UkSZK6w3jHo1VlPC41qDRJJ78pJEmS1PkmejxaVbbjUoNKE3T6m0KSJEmdbbLHo1VlOi41qDRYN7wpJEmS1LnqPR6tKstxqUGlgbrlTSFJkqTOde6KtXUfj1YNDiXnrljboIrq0/SgEhEfjYgNEZERsaiYNisiromIOyPi9oj4ekTMr1lnv4i4LiJ+FhE/iojnT2Reu3XLm0KSJEmd64KlC+nviylto78vuGDpwgZVVJ9WnFG5CjgG+OWw6ZcAh2TmM4EvAp+omfc+4JbM/GPg9cC/R8T0Ccxrq255U0iSJKlzLVkwl4tOW1z3cWl/X3DRaYtZsmBugyubnKYHlcz8VmZuGjbtkcz8amZWTz/cAhxUs8jLgYuLZb8H3AMcO4F5bdUtbwpJkiR1tnqPS8t0PFqWPip/SeWsChExG5iembUdNTYAB441r0V1jqsb3hSSJEnqfJM9Li3b8Wjbg0pEvAuYD5zb4O2eExGbqo+tW7c2cvNj6vQ3hSRJkrrDRI9Ly3g82tagEhFvB5YCL87M3wNk5n3AYETUvkoHARvHmjd825l5YWbOqz4GBgaa1YwRdfKbQpIkSd1jvOPSsh6Pti2oRMQ5wCuBF2bmb4fN/jywrFjuWcCTgG9OYF6pdOqbQpIkSd1ltOPSMh+PtuL2xMsjYhMwD1gZEesjYh7wQWAfYFVErImIW2tWewfw3Ij4GXAZ8OrM3D6BeaXTiW8KSZIkdZ/qcWlV2Y9H+5u9g8w8fZRZo14TlZn/BfyPyc4rq+qb4vRP3waU/00hSZKk7rRkwVz233sWW373SOmPR9vemb5XVN8U04LSvykkSZLUvQZm9vPUOQOlPx5t+hkVPWZgZj8DHfCmkCRJktrNMyqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCqSJEmSSsegIkmSJKl0DCpqqJXrNrP4vV9n5brN7S5FkiRJHcygooZZuW4zZ16xmvsf2saZV6w2rEiSJKluBhU1RDWkDA4lAINDaViRJElS3QwqmrLhIaXKsCJJkqR6GVQ0JaOFlCrDiiRJkuphUFHdxgspVYYVSZIkTVbTg0pEfDQiNkRERsSi8aYX8zZExB0RsaZ4nFoz748j4tsRcWdEfC8iFjS7DdrdRENKlWFFkiRJk9GKMypXAccAv5zg9KpTM3NR8biyZvpy4JLMPBh4P3BZg+vVOCYbUqoMK5IkSZqopgeVzPxWZm6a6PSxRMR+wJHAZ4pJVwNPjoj5U69UE1FvSKkyrEiS1NkcM02tMmZQiYhpEbHnCNP3jIhmh5zLI2JtRFwaEXOKaU8G7s3MQYDMTGAjcOAINZ4TEZuqj61btza53N5w7oq1dYeUqsGh5NwVaxtUkSRJahXHTFMrjRc2LgBePcL004p5zfL8zDwMWAz8GvjUZDeQmRdm5rzqY2BgoOFF9qILli6kvy+mtI3+vuCCpQsbVJEkSWoFx0xTq40XVP4EuHSE6Z8ETmp8ORWZubH4dzvwYeB5xaxfAftHRD9ARASVsykbm1WLdrVkwVwuOm1x3WGlvy+46LTFLFkwt8GVSZKkZnHMNLXDeEGlLzN3DJ9YTBtqRkER8biI2Kdm0iuBHxT73QKs5rGzPKcAmzJzfTNq0cjqDSuGFEmSOo9jpqldxgsqAxExc/jEYtpufVdGEhHLI2ITMA9YGRHrx5oOPBFYFRE/jIi1wLHAa2s2eTpwekTcCbwTeP1E6lBjTTasGFIkSeo8jpmmdhovqHwF+Gj1UiuA4ucPAddOZAeZeXrRT6Q/M5+YmfPHmX5XZh6emYdl5sLM/NPM3FCzvTsy8+jMPDgzj8xMe2W3yUTDiiFFkqTO45hparfxgsq7gUOAuyLiCxHxBeDnwNOBdzW7OJXfeGHFkCJJUudxzDSVwZhBJTMfyszjqFxe9e3i8frMPC4zH2pBfeoAo4UVQ4okSZ3HMdNUFhMaCyUzr8/MfyoeNzS7KHWealipMqRIktSZHDNNZdE/1syI+AEw6js1MxePNk+9Z8mCuey/9yy2/O4RQ4okSR3qgqULp3RGBRwzTY0xZlABzm5FEeoeAzP7GZgzYEiRJKlDVa+SqDeseFWFGmXMoJKZ32xVIZIkSSqHesOKIUWNNGYflYjYIyLeGhGviIi+iPjniFgbEVdFxJNaVaQkSZJayzHT1G7jdab/OHAS8BZgJbAP8P8CvwAubmplkiRJaivHTFM7jddHZXFmPj0iZgH3Aksycwj4WjFqvCRJkrrYeJeBGVLULOOdUXkUIDMfATYUIaVqe9OqkiRJUmk4ZpraYbygskdELIyIw4BZ1Z+L53u0oD5JkiSVgGOmqdXGu/RrJvBFoBqfv1Qzb2ojAUmSJKmjOGaaWmm8oPJp4HrgO5m5rQX1SJIkqcQcM02tMt6lX3OA5cBvIuIbEfH/RcRzImJaC2qTJEmS1KPGDCqZ+dbMPBSYD1wKHAhcTiW4XNuC+iRJkiT1oPEu/QIgM++NiKup3KJ4M/BKYFET65IkSZLUw8YMKhHxfOA44HjgScAtwLeAkzLzZ02vTpIkSVJPGu+Myo1UwsnfZeZ1zS9HkiRJksbvTP984KvA2yNifUT8e0S8JSIObkFtkiRJknrUmGdUMvNm4Gbg7yNiBnAUlcvAvhQRA5k5rwU1SpIkSeoxE+pMHxEHUAkoxwEnAPtRCTCSJEmS1HDjdab/OHAscADwHWAV8Grge5k52PzyJEmSJPWi8c6o/Ap4I3BLZm5vQT2SJEmSNO6Aj3+XmTcZUiRJktQLVq7bzOL3fp2V6za3u5SeN95dvyRJkqSesHLdZs68YjX3P7SNM69YbVhpM4OKJEmSel41pAwOJQCDQ2lYaTODiiRJknra8JBSZVhpL4OKJEmSetZoIaXKsNI+BhVJkiT1pPFCSpVhpT0MKpIkSeo5Ew0pVYaV1jOoSJIkqadMNqRUGVZay6AiSZKknlFvSKkyrLSOQUWSJEk949wVa+sOKVWDQ8m5K9Y2qCKNxqAiSZKknnHB0oX098WUttHfF1ywdGGDKtJoDCqSJEnqGUsWzOWi0xbXHVb6+4KLTlvMkgVzG1yZhjOoSJIkqafUG1YMKa1lUJEkSVLPmWxYMaS0nkFFkiRJPWmiYcWQ0h4GFUmSJPWs8cKKIaV9DCqSJEnqaaOFFUNKexlUJEmS1POqYaXKkNJ+BhVJkiSJSljZf+9ZTAsMKSXQ3+4CJEmSpLIYmNnPwJwBQ0oJeEZFkiRJUukYVCRJkiSVjkFFkiRJUukYVCRJkiSVjkFFkiRJUukYVCRJkiSVjkFFkiRJUukYVCRJkiSVjkFFkiRJUukYVCRJkiSVjkFFmoSV6zaz+L1fZ+W6ze0uRZIkqasZVKQJWrluM2desZr7H9rGmVesNqxIkiQ1kUFFmoBqSBkcSgAGh9KwIkmS1EQGFWkcw0NKlWFFkiSpeQwq0hhGCylVhhVJkqTmMKhIoxgvpFQZViRJkhrPoCKNYKIhpcqwIkmS1FgGFWmYyYaUKsOKJElS4xhUpBr1hpQqw4okSVJjGFSkGueuWFt3SKkaHErOXbG2QRVJkiT1JoOKVOOCpQvp74spbaO/L7hg6cIGVSRJktSbDCpSjSUL5nLRaYvrDiv9fcFFpy1myYK5Da5MkiSptzQ9qETERyNiQ0RkRCwab3ox748j4tsRcWdEfC8iFkxkntQI9YYVQ4okSVLjtOKMylXAMcAvJzgdYDlwSWYeDLwfuGyC86SGmGxYMaRIkiQ1VtODSmZ+KzM3TXR6ROwHHAl8pph0NfDkiJg/1rzmVK9eNtGwYkiRJElqvDL2UXkycG9mDgJkZgIbgQPHmSc13HhhxZAiSZLUHGUMKg0REedExKbqY+vWre0uSR1qtLBiSJEkSWqeMgaVXwH7R0Q/QEQElTMmG8eZt4vMvDAz51UfAwMDLWuAuk81rFQZUiRJkpqrdEElM7cAq4FXF5NOATZl5vqx5rW+UvWaJQvmsv/es5gWGFIkSZKarBW3J14eEZuAecDKiFg/1vTC6cDpEXEn8E7g9ROcJzXVwMx+njpnwJAiSZLUZP3N3kFmnj6Z6cW8O4CjJztPkiRJUnco3aVfkiRJkmRQkSRJklQ6BhVJkiRJpWNQkSRJapCV6zaz+L1fZ+W6ze0uRep4BhVJkqQGWLluM2desZr7H9rGmVesNqxIU2RQkSRJmqJqSBkcSgAGh9KwIk2RQUWSJGkKhoeUKsOKNDUGFUmSpDqNFlKqDCtS/QwqkiRJdRgvpFQZVqT6GFQkSVLLdMtdsSYaUqoMK9LkGVQkSVJLdMtdsSYbUqoMK9LkGFQkSVLTdctdseoNKVWd2m6pHQwqkiSpqbrprljnrlhbd0ipGhxKzl2xtkEVSd3LoCJJkpqm2+6KdcHShfT3xZS20d8XXLB0YYMqkrqXQUWSJDVFN94Va8mCuVx02uK6w0p/X3DRaYtZsmBugyuTuo9BRZKkEvGuWOVvd71hxZAiTY5BRZKkkvCuWN0bVgwp0uQZVCRJKgHvilXRSe2eaFgxpEj1MahIktRm3hVrV510V6zxwoohRaqfQUWSpDbyrli767S7Yo0WVgwp0tQYVCRJahPvirW7Tj24r7a7qlPbIZWJQUWSpDbwrli76/SD+yUL5rL/3rOYFnR0O6SyMKhIktRi3hVrd50eUqoGZvbz1DkDHd8OqQwMKpIktZB3xdpdt4QUSY1lUJEkqYW8K9auDCmSRmNQkSR1jG4Ytd27Yj3GkCJpLAYVSVJH6JZR270rVkWntkNS6xhUJEml1y2jtld5VyzviiVpfAYVSVKpddOo7bW8K1Znt0NS8xlUJEml1W2jtg/nXbEkaXQGFUlSKXXjqO0j8a5YkjQyg4okqXS6edT2kXhXLEnanUFFklQqvTBq+0i8K5Yk7cqgIkkqjV4atX0k3hVLkh5jUJEklUavjdo+Eu+KJUkVBhVJ6iKdPnJ7L47aLkkamUFFkrpEN4zc3qujtkuSdmdQkaQu0E0jt/fqqO2SpF0ZVCSpw3XjyO29Omq7JOkxBhVJ6mDdPHK7o7ZLUm8zqEhSh+qFkdsdtV2SepdBRZI6UC+N3O6o7ZLUmwwqktRhenHkdkdtl6TeY1CRpA7SyyO3O2q7JPUWg4okdZBeH7ndUdslqXcYVCT1FEdud+R2SVJnMKhI6hmO3G7fDklS5zCoSOoJjtxuSJEkdRaDiqSu58jthhRJUucxqEjqao7cbkiRJHUmg4qkruXI7YYUSVLnMqhI6kqO3G5IkSR1NoOKpK7jyO2GFElS5zOoSOoqjtzuyO2SpO5gUJG00/otW1n+zZ+zfsvWdpdSN0dud+R2SVJ3MKhIAmDVT7fw1bX38vD2HXx17b2s+umWdpdUF0dulySpOxhUJLFy3WbecfUPGcrKmYihTN5x9Q879vInR26XJKnzGVSkHlft07Fj2OVSOzq8r4Yjt0uS1NkMKlIPczDExxhSJEkqF4OK1KMcDPExhhRJksrHoCL1IAdDfIwhRZKkcjKoSD2mlwdDdOR2SZI6h0FF6iG9PhiiI7dLktQ5DCpSD+n1wRAduV2SpM5hUJF6iIMhOnK7JEmdwqAi9RAHQ5QkSZ2i6UElIj4aERsiIiNiUc30P46Ib0fEnRHxvYhYUDNvQ0TcERFrisepE1lP0vgcDFGSJHWCVpxRuQo4BvjlsOnLgUsy82Dg/cBlw+afmpmLiseVk1hP0jgcDFGSJJVd04NKZn4rMzfVTouI/YAjgc8Uk64GnhwR88faVr3rSdqdgyFKkqQya1cflScD92bmIEBmJrAROLBmmcsjYm1EXBoRcyaxHgARcU5EbKo+tm7d2sz2SB3JwRAlSVJZlbUz/fMz8zBgMfBr4FOT3UBmXpiZ86qPgYGBhhcpdYNqWJk2LKxMM6RIkqQ2aldQ+RWwf0T0A0REUDkrshEgM6v/bgc+DDxvIutJqs+SBXN5/ymH7XzeF8H7TznMkCJJktqmLUElM7cAq4FXF5NOATZl5vqIeFxE7FOz+CuBH4y3XksKl7rY8Yfux8DMfgI4ceH+HH/ofu0uSZIk9bD+Zu8gIpYDJwFzgZUR8WBmzgdOBy6LiHcBvwNeX6zyRODqiJgGBHAX8NqaTY62nqQpmtHfx4z+Gczfz0slJUlSezU9qGTm6aNMvwM4eoTpdwGHj7G9EdeTJEmS1D3K2plekiRJUg8zqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNIxqEiSJEkqHYOKJEmSpNKJzGx3DS0REY8C/93uOoABYGu7i2iybm+j7et83d5G29f5ur2Ntq/zdXsbu719UI42zsnMmaPN7JmgUhYRsSkz57W7jmbq9jbavs7X7W20fZ2v29to+zpft7ex29sHndFGL/2SJEmSVDoGFUmSJEmlY1BpvQvbXUALdHsbbV/n6/Y22r7O1+1ttH2dr9vb2O3tgw5oo31UJEmSJJWOZ1QkSZIklY5BRZIkSVLp9Le7AEmSJEmNEREbgEeBh2smLwTWAjOAQ4qfAe7IzFMj4jzg5cAOYCbw5cz865YVPQrPqNQhIjZExB0RcXtErI+IL0bEc4t5r4uIa2qW/URErCke24r1qs/3GmX7B0XEjRHxQESsaU2rdtl/s9t3QkR8NyJ+HBHrIuIDEdHS92IL2nh0zTLrImJ5RIw6oFGjNbt9NetGRNwQEb9tbot222+zf3/HRcTDNcutiYg9WtS8ag1N/x1GxMLib81PisfSFjStuu9m/w5fP+z39+uIWNGi5rWifX0RcWFU/o7+MCJWRcT8FjWvWkMr2vjPEfGjiPhpRFwaETM6vE1jfr5HxBsj4mcR8fOI+HhETJ9gzYtGmP6KiPhesb3vR8RNEXFKRCyrqfP+iLi75vnxY+wnh9dc/D/LiDi7ZtoREXFdRNxV7Pc/I+LkYt5lxfKH1yy/V0Rsrd12RAxExIeL38HtEfGD4r0w4uvRwtdgWfH/bU3xnryimP7ViPiLEZa/PSKWRuUzZc1uG2yuUzNzUc1jWmYuAk4EHqyZfmpE/C/gxcCzMvOZwDOAz7S43pFlpo9JPoANwKKa50uB3wJHAa8DrpnIemNsf1/gGOAkYE0Xtu9w4KnFz7OAm4HXdVkb9wSmFz/3AV8A3tYt7atZ/hzg48Bvu+z3d1w7/u+1uI17AncBxxTPp1EZIbgr2jfCej8CTumW9gEnA7fW/J15N/C5VrWvRW18M3ADlW+Ao/hb89cd3qZRP9+BpwD3AHOL9n4JOHOyNRfT3gT8FHh6zbRDhr9+wGXA2RN8bRL4PnBEzbSbgO9VtwEsAH4NvLRmmQOAP6/Z3/eBfxlW6/eqr0fR9huL3/cexbTpwDJgoF2vAXAklb+Z+9bUubj4+RTgthGW31LUftzw33cr38fD5h3EsM9s4GzgWoqbbJXp4RmVBsjMFcDFwNsbtL37M/Nm4KFGbG+qmtC+H2TmXcXPjwBrqPzHaZsmtPH3mbm9eDoD2IPKH/m2aHT7ACJiAZWDpfc1apv1akb7yqYJbXwVcEvxt4bM3JGZ/92gbU9aM3+HEXEUsB+VA7+2aEL7ksrlGbMiIoDHA5satO36Cmp8G58JfCMzt2XlaOprwGsatO0JafHn+/8CvpSZm4v2Xgy8ss5dnU/l4PvHNfu+IzP/qc7tVf0b8AaAiDiYykH4upr57wQ+mZlfrtnvPZn5qZplVgAviceuMng98Mma+ScA86mEtIeLbWzPzIszc+skaj2fxr4G84AHiwdZsbqY9yXgyRFxWM3ybwAurzkWaLUrY9ezymNdFfBZKkH5roi4PCLeMM7yLWNQaZxbqXyT0K2a0r6ImEvlj/NXGr3tOjS0jcUp/tupfLv0APCxRm27Tg1rX3H6/ePA6VSuZy2DRr9H/ygiVheXDZzRwO1ORSPb+HTg0Yj4SvEhdnlEzGnQtuvVrL+jbwQ+3cYDhqpGtu/LVL513gzcC/wJcF6Dtj0VjWzjbcD/jIjHF39zXk57vtRq1ef7gcAva55vKKZNSkTsBzyJSt2NtgI4MSJmUTkQ/7dh848AvjPONn4PfB04OSIOpXJm4ifDtnFbZm6rt8gmvQb/QSWkbIyIKyPiLyLiCVAJUsCneSzEzaISMi9t4P4na/ilXw+PtmBmbqbSh+U0Kn1XzgC+HU2+1HIiDCqNE+0uoMka3r6IeDyVD9sPZOb3G739OjS0jZm5ISvXes6l8s1ny67/H0Uj2/e3wIrM/Mm4S7ZOI9u3GpiXmYuBPwOWRcTLG7j9ejWyjf3AC6iEzcOBu4H/08Dt16MZf2ceB7yC9h4wVDWyfUdSuY78SVQurbmeyjfw7dbINl4GXAd8s3jcCQw2cPsT1fGf71Hpw7Q2Iu6Y4qYeBlYCLyse/3+d2/kklS8Q3sjuYacppvoaZObvgedR6ePxn1Q+038YEfsWi1wKnFYc3C8FflKyz8gxFWfVv12ccfp/qJxheUabyzKoNNCzqFwD3a0a2r6odDS8DvhiZpZlZNSm/A6LU9WfpfJNRTs1sn3HAv87KncWuRl4fNGZsZ3fyDesfZn5u8x8oPh5E5UP4+c1YttT1Mjf4UZgVWbeXVxm8hngOQ3adr2a8X/wZcC62ss/2qiR7XstcENm/jYzh4BPAaN2Am6hRv4/zMw8PzMPz8znAj9m18uMWqVVn+8bgT+seX5QMW1SMnMLlS8enl0z7XjgpcATp1YiUAkWFwLfzszfDZt3G3D0BGq8hUrAfgWVz8fh21g8lW/zm/UaFO/JH2TmR6mcxdxKpf8Jxd+Y9cU+3kA5vhyZkIg4MiL+qGbSoVQu6/tVm0rayaDSABHxp8BbgQ+2u5ZmaHT7ImKASki5LjP/vhHbnKomtHF+9e4kxR/bPwN+2Iht11lPQ9uXmc/LzD/MzIOodAz9XWYe1K4+Dk34/e0fxZ3oilD9EuAHjdj2FGpq9N+ZzwHPKs5sQuVbwtsbtO1Ja+Lf0TdSggOGJrTvLuCEmoO5l9DmL8ua8P9wVvXSmoj4Ayr9Hz7QiG1PooZWfr5fTeVSt7lFv6Nl7H4QP1F/B3youLSq6nFTLRAgM28F/h64YITZHwDeEBEnVScU7fnzEZb9S+DtmfngsOk3AL8APlpcQkVE9EfEW4rjh4lq6GsQEYcO64PyZGAOlf+LVZcC76ISkK6sd18NMryPylhfZMwGrijuZPYDKu14VTv7LVY5jkr9royIR6i86X8MnJiZt0bE04AlEVHbqfFzmXnORDccEXtSOcU9E9i72NanM/PcBtY/nqa1j8ofp2cDj4vHbof6+cz8h4ZUPnHNbOMJwFkRsYPK/7Prgfc2qvAJamb7yqCZ7TsFeGtEDFL5/X2eFl2eMEzT2piZGyPiH6lchzxE5dvHtzSy+Alo6ns0Ig4BFlEJYe3QzPZdBDwNuD0itlPpq7KsUYVPQjPbuDdwY/H+7AM+UttJu4na8vmemXdFxN9SuawIKn2Qlk9w0yuL90HVc6h02P9MROwN/DfwCHDmRGsdS2Z+ZJTpayPixcA/RMS/FDU8yAg3XcnM60fZRhZB5x+AdRHxMJXf/7VFG0bT7NdgTyrBZy6VS+ACeGdmrqlZ5krgw8CVk+z431DFl4ijzdsA7DNs2koql/SVTlTO+EuSJElSeXjplyRJkqTS8dKvNilunfcfI8z6emb+davrabRubx90fxttX+fr9jbavs7XjW3s5DZFxJfY/ZbIvyk6ovcEX4Ny8dIvSZIkSaXjpV+SJEmSSsegIkmSJKl0DCqSJEmSSsegIkklFBFPjYjri58/HxGL69zOpRHx44j4wgjzboyIX9QMCPbVKdS7KCJeUe/6jRIRX4mIo0aYnhGxZti01xfTz65jP1dFxOsmsNyNEXHysGkH1Lzm6yPi4ZrnH5psLZLUrbzrlySV0wuBb0TENOBwYM1kNxARTwReATw+M3eMstjbMvOaeoussQg4mTpG0o6I/swcnGoBxajVTwO+O8oigxFxRGbeVjx/A/D9qe53sjLzHiqvFxFxHPDhzFzU6jokqew8oyJJJRIRp0fELcA/Aq+iElCeQGUE+deMss5rIuKHxePaiHhSROwDrAJmAbdFxDsnUcOSiLg5Im6LiO9GxPHF9LkRsaqYvi4i/jUi+orbsf4dcHxxVuDiYvks6qhu99cRcVDx84aIeH9EfBf4VERMj4j3FftbExGfi4gnFMu+qTgrtCYi1o50xqTwYmBljn47y3+jEk6IiIOB6cC6mvoGIuKTEfGj4vG3NfMOjYhvF+2+Bnh8zby9IuLjRe0/jIhLImLGRF/vmu38a0S8q+b5IRHxq4joj4jzI+LqiLghIn4aEV+OiNnFcqO+dpLUyQwqklQimbk8M58D/IbKmZRPAO/NzOdk5qeHLx8RzwD+CXhxZh4GfBv4RGb+FjgReDAzF2Xm+0bZ5YdqLjs6MyKeCpwPnJiZR1AJS/8eETOB3wIvLaYfBhwEvDwztwDnAauKfS2bYHNnA0dl5mnAXwMPZeazi7MLa4G/L5b7IPAnxfTF1ISLYf4MuGaM/a0AToyIWVQCy78Nm/83wMyibUcBJ0fEqcW8TwOXZuaCYrlja9b7IHBTZj4beCaVz9a/HKOO0fwL8JbiLBrAGcAlNWebnge8KjMPBX4FXFBMH+u1k6SO5aVfklQyETEP2JKZgxFxBPDxMRY/HrguM+8unn8MOK/mYHc8u1z6FRFnAPOBb0VEdfIQlQHQ7gbeHxHHAAHsB/yIOi73KlxWc/bjZGDviDileD4D2FD8fD3w6Yj4MvC1zLxz+IYiYjrwXODPx9jfw8BK4GXF43Dg6Jr5LwD+KjOHgIci4nLghRHxNSqXal0GkJlrI+LmmvVOBo6OiHOK53sAo11qN6rMvCMifgz8aUSsBF4JLKxZ5NrM3Fz8fAmV4FXd/2ivnSR1LIOKJJVERDwZ+DKwN/C4ovP3wcDhEbE+M/9sApuZ6ii+QWUE7VeNUN+7qYSTozLzkYi4kMqlZaPZAdQGpuHLbh223/+dmSON6H0KcARwHPDViHh3Zg4PRycA/5mZ28eoBypnUb5CJdz9riaMjWSs17J2XgCnjBSg6vAR4B3AHCq/h/+aQA1jvXaS1LG89EuSSiIzf1VcuvM14NVULmW6JTMXjhFSVgEviogDiufLgOvH6Dw/npXACyLisOqEiHh28eMTgM1FSJlL5axE1e+oBKxa66lcQkVELAUeN8Z+rwHeFhF7FsvvGRELIqIf+KPM/H5m/jNwFfDsEdY/GdjtzmbDZeatVC6LumCE2d8A3hgVjwNeA/xHZv4O+AHw2qK2BcAxw2p/R1ErEfGEiJg/Xi2j+A9gLvBu4F+HzTsxKjdIAHhTUW91/7u9dnXuX5JKw6AiSeVzLHATxZ2/xlowM39EpY/CdRHxQyr9GN5c744zcz2VfinLI+L2iPgJcHYx+yPAURGxjkqfjdrargdmFp3JLy6mvQ34SESspnKZ1X1j7Pr9wPeAW4t23ELlcqtpQLWD+xoqZ1YurF0xKqdFlgDXTbCNH8nMH48w673Adip9PG4FvpSZnyvmvZZK/5EfUQk636pZ721ULitbU9R+PZX+O5NWXAp3KZVL/74zbPZNVPoL/RT4Q6Da8X60106SOlqMfnMUSZLKLyKeA7w7M1/S7loaISK+AlxZe/OEiDgf2Cczz25XXZLUap5RkSR1tMy8pRtCSkQcGRHrqdy84N/bXY8ktZtnVCRJkiSVjmdUJEmSJJWOQUWSJElS6RhUJEmSJJWOQUWSJElS6RhUJEmSJJWOQUWSJElS6fxfDQucliPTOSgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 960x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = [\"DLT_{}\".format(x) for x in range(1, NUM_OF_REGRESSORS + 1)] + ['LGT_MCMC', 'LGT_SVI','ETS']\n",
    "fig, ax = plt.subplots(1, 1,figsize=(12, 6), dpi=80)\n",
    "markerline, stemlines, baseline = ax.stem(\n",
    "    np.arange(len(labels)), np.array(WBIC_ls), label='WBIC', markerfmt='D')\n",
    "baseline.set_color('none')\n",
    "markerline.set_markersize(12)\n",
    "ax.set_ylim(1020, )\n",
    "ax.set_xticks(np.arange(len(labels)))\n",
    "ax.set_xticklabels(labels)\n",
    "# because list type is mixed index from 1;\n",
    "ax.axvline(x=NUM_OF_EFFECTIVE_REGRESSORS - 1, color='red', linewidth=3, alpha=0.5, linestyle='-', label='truth') \n",
    "ax.set_ylabel(\"WBIC\")\n",
    "ax.set_xlabel(\"# of Features / Model Type\")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BIC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we use DLT model as an example. Different DLT models (the number of features used changes) are fitted and their BIC values are calculated respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-03-25T01:06:52.190878Z",
     "start_time": "2022-03-25T01:06:51.192123Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n",
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n",
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n",
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BIC value with 1 regressors: 1247.445\n",
      "------------------------------------------------------------------\n",
      "BIC value with 2 regressors: 1191.889\n",
      "------------------------------------------------------------------\n",
      "BIC value with 3 regressors: 1139.408\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n",
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BIC value with 4 regressors: 1079.639\n",
      "------------------------------------------------------------------\n",
      "BIC value with 5 regressors: 1082.445\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n",
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BIC value with 6 regressors: 1082.482\n",
      "------------------------------------------------------------------\n",
      "BIC value with 7 regressors: 1082.217\n",
      "------------------------------------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n",
      "INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BIC value with 8 regressors: 1081.698\n",
      "------------------------------------------------------------------\n",
      "BIC value with 9 regressors: 1078.876\n",
      "------------------------------------------------------------------\n",
      "BIC value with 10 regressors: 1187.269\n",
      "------------------------------------------------------------------\n",
      "CPU times: user 935 ms, sys: 58.7 ms, total: 994 ms\n",
      "Wall time: 995 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "BIC_ls = []\n",
    "for k in range(0, NUM_OF_REGRESSORS):\n",
    "    regressor_col = x_cols[:k + 1]\n",
    "    dlt_mod = DLT(\n",
    "        estimator='stan-map',\n",
    "        response_col=response_col,\n",
    "        date_col=dt_col,\n",
    "        regressor_col=regressor_col,\n",
    "        seed=2022,\n",
    "        # fixing the smoothing parameters to learn regression coefficients more effectively\n",
    "        level_sm_input=0.01,\n",
    "        slope_sm_input=0.01,\n",
    "        )\n",
    "    dlt_mod.fit(df=df)\n",
    "    BIC_temp = dlt_mod.get_bic() \n",
    "    print(\"BIC value with {:d} regressors: {:.3f}\".format(k + 1, BIC_temp))\n",
    "    print('------------------------------------------------------------------')\n",
    "    BIC_ls.append(BIC_temp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot below shows the BIC vs the number of features (blue line). The true model is indicated by the vertical red line. The horizontal gray line shows the minimum (optimal) value. The minimum is at the true value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-03-25T01:06:52.363849Z",
     "start_time": "2022-03-25T01:06:52.193136Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyoAAAGbCAYAAAArwGQ+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAxOAAAMTgF/d4wjAAAnvUlEQVR4nO3df5TdZX0v+vcnmQBiEDSFBhsp9UTURjAEW69d2FatJ63acyxcRUV7dfV64MgqVWqvxevtaY/3HrTrFH/c4iVaPWjl3ENFjrbVEq1I1esSlQjGVMVoYwySUqFigwch5Ll/zN4yhMzMnsnes78z83qttVdmP8/+Pvv5ftZk9vf9/bWrtRYAAIAuWTHuCQAAABxMUAEAADpHUAEAADpHUAEAADpHUAEAADpHUAEAADpnYtwTWChHHnlkO/7448c9DWCcDhxI9u174Pnq1ckK+2sAYBxuvfXWe1trR07Xv2yCyvHHH589e/aMexrAON1+e/KOdzzw/FWvSk44YXzzAYBlrKr+aaZ+uxIBAIDOEVQAAIDOEVQAAIDOWTbXqAAAwEI7cOBAWmvjnsbYVFVWzPPGNYIKAAAM2b333pvdu3fnvvvuG/dUxm7VqlU56aSTcsQRR8xpOUEFAACGbPfu3TnmmGOyZs2aVNW4pzM2rbXccccd2b17d9avXz+nZQUVAAAYogMHDuS+++7LmjVrMjFhc3vNmjW58847c+DAgTmdBjbyi+mr6u1VtauqWlVtnNL+sar6clXdVFWfrqrTp/Q9rqo+W1W3VNUXqmrDIH0AADBu/WtSlvORlKn6dZjrtToLcdevq5OcmeTbB7W/sLV2WmttY5JLk1wxpW9Lkne21k5J8uY59AEAAEvAyINKa+1TrbWHfCV8a+37U54em6QlSVWdkOQpSd7f6/tgksdU1fqZ+kY0fQAAWBL+8A//MPfcc8+8ln3rW9+avXv3PmisV7/61UOa2aGN9XtUqup9VfWdJG9M8rJe82OS3NZa258kbfIY0e4kJ83SBwAATOOP/uiPDhlU9u/fP+uyBweVhTDWq3taa7+ZJFX1v2TyNK7nDGvsqrooyUX958cee+ywhgYAgMHt35/ceefo3+dRj0qmuXj//PPPT5I8/elPz8qVK/PoRz86a9euzc6dO3P77bfna1/7Wqoq//zP/5zjjjsuSfITP/ET+eIXv5j3ve99+e53v5tzzjknD3vYw3LFFVckSW677bb8+q//er75zW9m7dq1ufrqq/OoRz1qaKvTidsQtNbeW1WXV9WaJN9JcmJVTbTW9tfk1TcnZfLIyQ9m6Dt4zEszee1LkmTdunXL95t2AAAYnzvvTN7xjtG/z6telZxwwiG7Lr/88mzZsiWf/vSnc9xxx+XlL395brzxxnzmM5/JMcccM+Owf/AHf5D3vOc9ueqqq7Jx48YkyYc+9KHccMMNufHGG7NmzZq86EUvypYtW3LxxRcPbXXGcupXVR1XVY+e8vz5Se5Icmdr7fYk25K8tNd9dpI9rbWdM/Ut2OQBAGAJeMELXjBrSJnJr/7qr2bNmjVJkqc97Wn55je/OaypJVmAIypVtSXJc5OsTbK1qv4lybOSfKCqHpbkQJJ/SvK89sA9y85LckVVvT6TR1FeMWXImfoAAIABrF69+kHPV65cmfvvv//Hz2e78P6oo4560LKDXOsyFyMPKq2186bp+vkZlvl6kqfNtQ8AADrnUY+aPC1rId5nBsccc0zuuuuuH1+DcrD169fnhhtuyHOe85xcc801ufvuu3/c94hHPCJ33XXXMGc7q05cowIAAEvWxMS0144spN/93d/Ns5/97Bx99NF59KMf/ZD+t7zlLbnwwgvzhje8Ic997nN/fFpXklx44YV55StfmaOPPvrHF9OPWs31GyIXq3Xr1rU9ex7ydS7AcnL77Q++mHGGiw4BYL7uv//+3HLLLTnllFOycuXKcU9n7KarR1Xd2lpbN91yY/0eFQAA4MG27tibTW/8eLbuWNjvLekaQQUAADpi6469ueDKbbnz7ntzwZXblnVYEVQAAKAD+iFl/4HJSzP2H2jLOqwIKgAAMGYHh5S+5RxWBBUAABij6UJK33INK4IKAACMyWwhpW85hhVBBQAAxmDQkNK33MKKoAIAAAtsriGl73DCysknn5zHP/7x2bhxY574xCfmJS95Se6+++5cf/312bhx449ft2/fvrz61a/O+vXr8+QnPzmnn356Xvva1+a+++6b83seDkEFAAAW0HxDSt/hhJWrrroqN910U3bs2JG77rrrId8y31rL8573vNx9993Zvn17br755nz+85/P+vXr86Mf/Whe852viQV9NwAAWOYuvmb7vENK3/4DLRdfsz2bN6yd1/L33ntvfvjDH+aRj3zkg9qvu+667Ny5Mx/72MdyxBFHJElWrVqV888//7DmOx+OqAAAwAK65KxTM7GiDmuMiRWVS846dc7LnXPOOdm4cWPWrl2bFStW5IUvfOGD+m+88cacccYZPw4p4ySoAADAAtq8YW0uO3fTvMPKxIrKZedumtfRlP6pX9/73vdy8skn53Wve9285rAQBBUAAFhg8w0rhxNSHjTOxETOPvvsXHvttQ9qP+OMM7Jt27bce++9hzX+MAgqAAAwBnMNK8MKKX3XXXddHv/4xz+o7ZnPfGZ+5md+JhdeeGHuueeeJMn+/fvzzne+M/v27RvK+w5KUAEAgDEZNKwMK6T0r1F50pOelK9+9at529ve9qD+qspHPvKRHHHEEdmwYUOe9KQn5bTTTss3vvGNHHXUUYf13nNVrR3eHQcWi3Xr1rU9e/aMexrAON1+e/KOdzzw/FWvSk44YXzzAWBJuv/++3PLLbfklFNOycqVKwdaZqZbFg/7SMpCm64eVXVra23ddMs5ogIAAGM23ZGVxR5SDoegAgAAHXBwWFnOISURVAAAYKiqJoPGfC6x6IeVRz38iCUTUvp16NdlUL6ZHgAAhmjFihVZtWpV7rjjjqxZs2bOG+i/8oTj8yuvf2aSyes7FrPWWu64446sWrUqK1bM7RiJoAIAAEN20kknZffu3bnzzjvHPZWxW7VqVU466aQ5LyeoAADAkB1xxBFZv359Dhw4MK9TwJaKqprzkZQ+QQUAAEZkvhvpuJgeAADoIEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEEFAADoHEFlAW3dsTeb3vjxbN2xd9xTAQCAThNUFsjWHXtzwZXbcufd9+aCK7cJKwAAMANBZQH0Q8r+Ay1Jsv9AE1YAAGAGIw8qVfX2qtpVVa2qNvbajqqqD1XVLVV1c1V9vKrWT1nm+qr6h6q6qfd4zZS+E6rq2qr6RlV9pap+cdTrcDgODil9wgoAAExvIY6oXJ3kzCTfPqj9nUke31p7cpIPJ/mzg/pf01rb2Hu8ZUr7m5J8rrX2uCSvSPJfq2rViOZ+WKYLKX3CCgAAHNrIg0pr7VOttT0Htd3TWvtoa62/Bf+5JCcPOOQLk1zeG+cLSb6b5JeGNN2hmS2k9AkrAADwUF25RuV3MnlUZao3VdX2qrqqqh6bJFW1Jsmq1trUrfpdSU46eMCquqiq9vQf+/btG9XcH2LQkNInrAAAwIONPahU1euTrE9y8ZTml7XWnpDktCSfTvLXcx23tXZpa21d/7F69erhTHgWcw0pfcIKAAA8YKxBpapem+SsJL/WWvthv7219p3ev6219qdJHltVa1prdyTZX1VrpwxzcpLdCzjtac03pPQJKwAAMGlsQaWqLkry4iTPbq19f0r7RFX95JTnZyf5x15ISZIPJDm/1/dzSX4qyd8t1LxncvE12+cdUvr2H2i5+JrtQ5oRAAAsTgtxe+ItVbUnybokW6tqZ1WtS/InSY5L8sneLYhv6C1yZJKP9K5PuTnJq5L8mylDvi7JL1TVN5JckeSlrbX7Rr0eg7jkrFMzsaIOa4yJFZVLzjp1SDMCAIDFaWLUb9BaO2+arkNu0bfW7k7ylBnG+8ck/3oIUxu6zRvW5rJzN8379K+JFZXLzt2UzRvWzv5iAABYwsZ+Mf1S0w8rcz2yIqQAAMADBJURmGtYEVIAAODBBJURGTSsCCkAAPBQgsoIzRZWhBQAADg0QWXEpgsrQgoAAExPUFkA/bDSJ6QAAMDMBJUFsnnD2px47FFZWRFSAABgFiP/HhUesPrIiaw+frWQAgAAs3BEBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6BxBBQAA6JyRB5WqentV7aqqVlUbe21HVdWHquqWqrq5qj5eVeunLHNCVV1bVd+oqq9U1S8O0gcAACwNC3FE5eokZyb59kHt70zy+Nbak5N8OMmfTel7U5LPtdYel+QVSf5rVa0aoA8AAFgCRh5UWmufaq3tOajtntbaR1trrdf0uSQnT3nJC5Nc3nvtF5J8N8kvDdAHAAAsAV25RuV3MnlUJVW1Jsmq1treKf27kpw0U98CzRMAAFgAE+OeQFW9Psn6JM8a8rgXJbmo//zYY48d5vAAAMAIjfWISlW9NslZSX6ttfbDJGmt3ZFkf1WtnfLSk5Psnqnv4LFba5e21tb1H6tXrx7VajACW3fszaY3fjxbd+yd/cUAACw5YwsqvSMeL07y7Nba9w/q/kCS83uv+7kkP5Xk7wboYwnYumNvLrhyW+68+95ccOU2YQUAYBlaiNsTb6mqPUnWJdlaVTural2SP0lyXJJPVtVNVXXDlMVel+QXquobSa5I8tLW2n0D9LHI9UPK/gOT91nYf6AJKwAAy9DIr1FprZ03TVfNsMw/JvnXc+1jcTs4pPT1w8pl527K5g1rp1kaAIClpCt3/WKZmy6k9DmyAgCwvAgqjN1sIaVPWAEAWD4EFcZq0JDSJ6wAAByexXJ3VUGFsZlrSOkTVgAA5mcx3V1VUGEs5htS+oQVAIC5WWx3VxVUGIuLr9k+75DSt/9Ay8XXbB/SjAAAlq7Z7q7axbAiqDAWl5x1aiZWTHuH6oFMrKhcctapQ5oRAMDStFjvriqoMBabN6zNZedumndYmVhRvlcFAGAWi/nuqoIKYzPfsCKkAADMbrHfXVVQYazmGlaEFACA2S2Fu6sKKozdoGFFSAEAmN1SubuqoEInzBZWhBQAgMEslburCip0xnRhRUgBABjcUrm7qqBCp/TDSp+QAgAwN0vl7qqCCp2zecPanHjsUVlZ6cR/EgCAxWYp3F1VUKGTVh85kccev7oT/0kAABajxX53VUEFAACWqMV8d1VBBQAAlrDFendVQQUAAJa4xXh3VUEFAACWgcV2d1VBBQAAlonFdHfViXFPAAAAWDirj5zI6kVwd1VHVAAAgM4RVAAAgM4RVAAAgM4RVAAAgM4RVAAAgM4RVAAAgM4RVAAAgM4RVAAAgM4RVAAAgM4RVAAAgM4RVAAAgM6ZNahU1cqqOvoQ7UdXlaADAAAM3SBB45IkLz1E+7m9PgAAgKEaJKg8K8m7D9H+niTPHe50AAAABgsqK1pr9x/c2Gs7MPwpAQAAy90gQWV1VR15cGOv7SHXrgAAAByuQYLKXyd5e1VN9Bt6P78lyUdGNTEAAGD5mpj9JXlDJgPJt6rqxl7bpiT/ENeoAAAAIzBrUGmt3Z3kl6vqWZkMKEnyf7fWrhvpzAAAgGVrkCMqSZLW2ieSfGKEcwEAAEgyQFCpqi8laYfqStJaa5sO0QcAADBvgxxRefWoJwEAADDVIHf9+nySn01yYpJPJ/n1JH+a5LeT7Bzd1AAAgOVqkKDyrkze3evfJdma5Lgk/1sm7/p1+chmBgAALFuDnPq1qbX2s1V1VJLbkmxurR1I8jdVtX200wMAAJajQY6o/ChJWmv3JNnVCyl9941kVgAAwLI2yBGVo6rq1Eze5Wvqz0nysJHNDAAAWLYGCSoPS/KXU55P/flQty0GAAA4LIN8M/3JCzAPAACAHxvkGhUAAIAFJagAAACdI6gAAACdM/KgUlVvr6pdVdWqauNs7b2+XVX19aq6qfc4Z0rf46rqs1V1S1V9oao2jHodAACAhbUQR1SuTnJmkm8P2N53TmttY+9x1ZT2LUne2Vo7Jcmbk1wx5PkCAABjNvKg0lr7VGttz6DtM6mqE5I8Jcn7e00fTPKYqlp/+DMFAAC6osvXqLyvqrZX1bur6vhe22OS3NZa258krbWWZHeSkw5euKouqqo9/ce+ffsWbuYAAMBh6WpQ+cXW2mlJNiX5XpL3znWA1tqlrbV1/cfq1auHPkkAAGA0Bvlm+gXXWtvd+/e+qnprklt6Xd9JcmJVTbTW9ldVZfJoyu7xzBQAABiFzh1RqaqHV9VxU5penORLSdJauz3JtiQv7fWdnWRPa23ngk4SAAAYqZEfUamqLUmem2Rtkq1V9S+ttfXTtSf5ySQfrKqVSSrJt5L85pQhz0tyRVW9PskPkrxi1OsAAAAsrJEHldbaeXNs/1aS02cY7+tJnjac2cHit3XH3lx8zfZcctap2bxh7binAwAwFJ079QsY3NYde3PBldty59335oIrt2Xrjr3jnhIAwFAIKrBI9UPK/gMtSbL/QBNWAIAlQ1CBRejgkNInrAAAS4WgAovMdCGlT1gBAJYCQQUWkdlCSp+wAgAsdoIKLBKDhpQ+YQUAWMwEFVgE5hpS+oQVAGCxElSg4+YbUvqEFQBgMRJUoOMuvmb7vENK3/4DLRdfs31IMwIAGD1BBTrukrNOzcSKOqwxJlZULjnr1CHNCABg9AQV6LjNG9bmsnM3zTusTKyoXHbupmzesHbIMwMAGB1BBRaB+YYVIQUAWKwEFVgk5hpWhBQAYDETVGARGTSsCCkAwGInqMAiM1tYEVIAgKVAUIFFaLqwIqQAAEuFoAKLVD+s9AkpAMBSIqjAIrZ5w9qceOxRWVkRUgCAJWVi3BMADs/qIyey+vjVQgoAsKQ4ogIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHSOoAIAAHTOyINKVb29qnZVVauqjbO19/oeV1WfrapbquoLVbVhkD4AAGBpWIgjKlcnOTPJtwdsT5ItSd7ZWjslyZuTXDFgHwAAsASMPKi01j7VWtszaHtVnZDkKUne32v6YJLHVNX6mfpGM3sAAGAcuniNymOS3NZa258krbWWZHeSk2bpAwAAloguBpWhqKqLqmpP/7Fv375xTwkAABhQF4PKd5KcWFUTSVJVlckjJrtn6XuQ1tqlrbV1/cfq1asXbAUAAIDD07mg0lq7Pcm2JC/tNZ2dZE9rbedMfQs/UwAAYFQW4vbEW6pqT5J1SbZW1c6Z2nvOS3JeVd2S5PeTvGLAPgAAYAmYGPUbtNbOm0t7r+/rSZ421z4AAGBp6NypXwAAAIIKAADQOYIKAADQOYIKAADQOYIKAADQOYIKAADQOYIKAADQOYIKAADQOYIKAIe0dcfebHrjx7N1x95xT6Xz1GpwagUMSlAB4CG27tibC67cljvvvjcXXLnNRuUM1GpwajU3Qh3LnaACwIP0Nyb3H2hJkv0Hmo3KaajV4NRqboQ6EFSAZWbn7fuy5e++mZ237xv3VDrp4I3JPhuVD6VWg1OruRHqYJKgAiwbn/za7fno9tvyP+67Px/dfls++bXbxz2lTpluY7LPxtID1GpwajU3Qt38OE1uaRJUgGVh6469ed0Hv5wDbfLD/0Bred0Hv+xDrWe2jck+G0tqNRdqNTdC3fw4TW7pElSAJa//IXb/QR/+9/vQTzL4xmTfct5YUqvBqdXcCHXz4zS5pU1QAZY0eyhnNteNyb7lWDe1GpxazY1QNz9Ok1v6BBVgybKHcmbz3ZjsW051U6vBqdXcCHXzYyfU8iCoAEuSPZSzu/ia7fPemOzbf6Dl4mu2D2lG3aVWg1OrwQl182Mn1PIhqABLjj2Ug7nkrFMzsaIOa4yJFZVLzjp1SDPqLrUanFoNTqibOzuhlhdBBVhS7KEc3OYNa3PZuZvmvVE5saJy2bmbsnnD2iHPrHvUanBqNTihbm7shFp+BBVgSbGHcm7mu1G5nDYm+9RqcGo1GKFucHZCLU+CCrCk2EM5d3PdWFpOG0cHU6vBqdVghLrB2Am1PAkqwJJiD+X8DFq35VqfqdRqcGo1GKFudnZCLU+CCrDk2EM5P7PVbbnXZyq1GpxaDUaom5mdUMuToAIsSfZQzs90dVOfh1KrwanVYIS6mdkJtfwIKsCSZQ/l/PTr1qc+01OrwanVYIS6mdkJtbwIKsCSZg/l/GzesDYnHntUVlbUZxZqNTi1GoxQNzM7oZYPQQVY8vofaisP+lBb6UNsRquPnMhjj1+tPgNQq8Gp1WCEupnZCbU8CCrAsrB5w9q8+ezTfvx8RVXefPZpPsSAzhLqZuY0uaVPUAGWjWc84YSsPnIileQ5p56YZzzhhHFPCYDD4DS5pU1QAZaVIyZW5JEPPyLrT1g97qkAMAROk1u6JsY9AQAAOByrj5zIaqfJLTmOqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0jqAAAAJ0z8qBSVW+vql1V1apq45T2x1XVZ6vqlqr6QlVtmNK3q6q+XlU39R7nDLIcAACwNCzEEZWrk5yZ5NsHtW9J8s7W2ilJ3pzkioP6z2mtbew9rprDcgAAwCI38qDSWvtUa23P1LaqOiHJU5K8v9f0wSSPqar1M4013+UAAIDFZVzXqDwmyW2ttf1J0lprSXYnOWnKa95XVdur6t1VdfwclkuSVNVFVbWn/9i3b98o1wcAABiirl5M/4uttdOSbEryvSTvnesArbVLW2vr+o/Vq1cPfZIAAMBoTIzpfb+T5MSqmmit7a+qyuRRkd1J0lrr/3tfVb01yS2DLAcAACwNYzmi0lq7Pcm2JC/tNZ2dZE9rbWdVPbyqjpvy8hcn+dJsyy3IxAEAgAUx8iMqVbUlyXOTrE2ytar+pbW2Psl5Sa6oqtcn+UGSV/QW+ckkH6yqlUkqybeS/OaUIadbDgAAWCJGHlRaa+dN0/71JE87RPu3kpw+w3iHXA4AAFg6unoxPQAAsIwJKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOcIKgAAQOdUa23cc1gQVfWjJP807nkkWZ1k37gnsUio1eDUam7Ua3BqNTi1GpxaDU6tBqdWc9OFeh3fWjtyus5lE1S6oqr2tNbWjXsei4FaDU6t5ka9BqdWg1OrwanV4NRqcGo1N4uhXk79AgAAOkdQAQAAOkdQWXiXjnsCi4haDU6t5ka9BqdWg1OrwanV4NRqcGo1N52vl2tUAACAznFEBQAA6BxBBQAA6BxBBQAA6BxBZR6qaldVfb2qbq6qnVX14ar6hV7fy6vqQ1Ne+2dVdVPvcW9vuf7zY6YZ/+Squr6q7qqqmxZmrUZjAWr1zKr6fFX9fVXtqKo/rqpF+Xu9ALV62pTX7KiqLVU17ZcsddmoazVl2aqq66rq+6Ndo9FagN+tX66q/zHldTdV1cMWaPWGaiF+t6rq1N7f+K/2HmctwKoN3QL8Xr3ioN+p71XVNQu0ekO1ALVaUVWX1uRn4Zer6pNVtX6BVm+oFqhW/7mqvlJVX6uqd1fVEQu0eodtAeoz4zZoVf1WVX2jqr5ZVe+qqlWjWtckSWvNY46PJLuSbJzy/Kwk30/y1CQvT/KhQZabYfxHJTkzyXOT3DTu9e14rU5P8tjez0cl+UySl497vTtaq6OTrOr9vCLJf0/ymnGvdxdrNeX1FyV5V5Lvj3udu1yvJL+82P9WLWCtjk7yrSRn9p6vzOQ3M4993btWq0Ms95UkZ497vbtYqyTPT3LDlL/xb0jyF+Ne747W6pVJrktyRJLq/Y3/vXGvd4fqM+02aJKfSfLdJGt7tfvLJBeMcn0X5Z7nrmmtXZPk8iSvHdJ4d7bWPpPk7mGM1yUjqNWXWmvf6v18T5Kbkpw8jLHHbQS1+mFr7b7e0yOSPCzJkrjt37BrlSRVtSGTH/5vGtaYXTGKei1VI6jVS5J8rvc3Pq21+1tr/zSkscdqlL9XVfXUJCdkcsNo0RtBrVqSI5McVVWV5BFJ9gxp7LEaQa2enORvW2v3tsmt779J8rIhjb3gFngb9H9O8pettb292l2e5MXDeN/pCCrDc0OSDeOexCIxklpV1dpM/if662GPPUZDrVXvkO7NSb6X5K4k7xjW2B0wtFr1DmW/K8l5Se4fxpgdNOz/h/+qqrZV1Req6lVDHLcLhlmrn03yo6r6697pF++rquOHNHYXjOqz8LeS/PmUnS1LwTBr9VdJrk+yN8ltSZ6V5A+GNHYXDLNWNyb5N1X1iN7f+hdm8e/gXKht0JOSfHvK8129tpERVIanxj2BRWTotaqqR2TyD/Uft9a+OOzxx2iotWqt7WqtPTmTh22PzOQh46VimLX6D0muaa19dYhjds0w67UtybrW2qYkv5Hk/Kp64RDHH7dh1moiya9kMgSfnuTWJP/PEMcft1H8fX94khclefewxx6zYdbqKUmelOSnkjw6yScyubd7qRhmra5Icm2Sv+s9bkmyf4jjj8OS3QYVVIbn5zJ5/iyzG2qteheEXZvkw621zn/L6hyN5PeqtbYvyX9Lcu6wxx6jYdbql5L8dlXtyuR1T4/oXcC4lPZ8D61erbUftNbu6v28J8n/m+Tpwxi7I4b5u7U7ySdba7f2Tp14f5L/aUhjd8Eo/ma9IMmO1trfD3nccRtmrX4zyXWtte+31g4keW+SZwxp7C4Y5t+r1lr7w9ba6a21X0jy90l2DGPsMVqobdDdSX56yvOTe20jI6gMQVX92yT/PsmfjHsuXTfsWlXV6kyGlGtba//nMMbsihHUan3/7hy9O5z8RpIvD2PscRt2rVprT2+t/XRr7eRMXlT4g9bayUvlWoIR/G6dWL277fV2HDwvyZeGMfa4jeDv+18k+bneUeAkeU6Sm4c09liN8LPwt7LEjqaMoFbfSvLMKXevel6WyM7TEfy9OqqqHtn7+SeS/H6SPx7G2OOwwNugH8zkaXNre9dCnZ/JnZ4jMzHKwZe4q6rqniQPz2Qaf05r7YaqemKSzVU19SK2v2itXTTowFV1dCYPRR6Z5NjeWH/eWrt4iPNfSCOrVZLfSfLzSR5eD9zi8wOttf9rKDNfeKOs1TOTXFhV92fy//4nkrxxWBMfg1HWaikaZb3OTvLvq2p/Jn+3PpDkvwxr4mMwslq11nZX1X9K8tmqOpDJU7/+3TAnv8BG+v+wqh6fZGMmA91iN8paXZbkiUlurqr7MnmtyvnDmvgYjLJWxya5vvf/b0WSt7XW/mpoM18YY9kGba19q6r+Q5L/r/fy65NsGcL6TD+f3u3GAAAAOsOpXwAAQOc49WtMquqEJB87RNfHW2u/t9Dz6TK1GpxaDU6t5ka9BqdWg1OrwanV4NRqZoupPk79AgAAOsepXwAAQOcIKgAAQOcIKgAAQOcIKgAMpKoeW1Wf6P38garaNM9x3l1Vf19V//0QfddX1T9U1U29x0cPY74bq+pF810egPFy1y8ABvXsJH9bVSuTnJ7kprkOUFU/meRFSR7RWrt/mpe9prX2oflOcoqNSZ6feXxzclVNtNb2D2EOAMyTIyoAzKiqzquqzyX5T0leksmA8shMfrv6y6ZZ5mVV9eXe4yNV9VNVdVySTyY5KsmNVfX7c5jD5qr6TFXdWFWfr6pn9NrXVtUne+07qupPq2pF7/ab/zHJM3pHZi7vvb715tEf93tVdXLv511V9eaq+nyS91bVqqp6U+/9bqqqv6iqR/Ze+7/2jgrdVFXbq+qpcyoqALNye2IABlJVO5M8IckFmfz8eOs0r3tSkr9NckZr7daq+t+TnNla+7VeKLiptXbcNMten+Snk9zVa3pXkr9JcmWSza21H1TV+iSfTnJykkoy0Vrb1zvS8+Ek72+t/beqenmS57fWnj9l/Jbkka217/eefy/JU1pru6pqV2/er2yttap6fZKVrbU39l77fyRZ21q7oKruSvKE1tptVbUqyZGttX1zKCcAs3DqFwCzqqp1SW5vre2vqjMyGSCm84wk17bWbu09f0eSP+gFiUE86NSvqnpVkvVJPlVV/eYDSU5KcmuSN1fVmZkMLSck+UrmcbpXzxXtgT14z09ybFWd3Xt+RJJdvZ8/keTPq+qvkvxNa+2Web4fANMQVACYVlU9JslfJTk2ycOr6qYkpyQ5vap2ttZ+Y4BhDvfQfWXyG5Nfcoj5vSGT4eSprbV7qurSTJ5aNp37k0wNTAe/dupRkUry2621Q32D89lJzkjyy0k+WlVvaK3NNxwBcAiuUQFgWq2177TWNmby9KuXJvmNJJ9rrZ06Q0j5ZJJfrapH956fn+QTM1w8P5utSX6lqk7rN1TVz/d+fGSSvb2QsjbJC6Ys94NMBqypdiZ5am+Ms5I8fIb3/VCS11TV0b3XH11VG6pqIsm/aq19sbX2n5NcneTnZxgHgHkQVAAYxC9l8rqQZ2fyOo5ptda+kuT3klxbVV9O8vQkr5zvG7fWdmbyIv4tVXVzVX01yat73W9L8tSq2pHkzw+a2yeSHNm7oP/yXttrkrytqrZl8s5ld8zw1m9O8oUkN/TW43OZvJPYyiTvqaqv9I4wnZHk0vmuHwCH5mJ6AACgcxxRAQAAOkdQAQAAOkdQAQAAOkdQAQAAOkdQAQAAOkdQAQAAOkdQAQAAOkdQAQAAOuf/B2vzW81spbMiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 960x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = [\"DLT_{}\".format(x) for x in range(1, NUM_OF_REGRESSORS + 1)]\n",
    "fig, ax = plt.subplots(1, 1,figsize=(12, 6), dpi=80)\n",
    "markerline, stemlines, baseline = ax.stem(\n",
    "    np.arange(len(labels)), np.array(BIC_ls), label='BIC', markerfmt='D')\n",
    "baseline.set_color('none')\n",
    "markerline.set_markersize(12)\n",
    "ax.set_ylim(1020, )\n",
    "ax.set_xticks(np.arange(len(labels)))\n",
    "ax.set_xticklabels(labels)\n",
    "# because list type is mixed index from 1;\n",
    "ax.axvline(x=NUM_OF_EFFECTIVE_REGRESSORS - 1, color='red', linewidth=3, alpha=0.5, linestyle='-', label='truth') \n",
    "ax.set_ylabel(\"BIC\")\n",
    "ax.set_xlabel(\"# of Features\")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. Watanabe Sumio (2010). \"Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory\". Journal of Machine Learning Research. 11: 3571–3594.\n",
    "2. McElreath Richard (2015). \"Statistical Rethinking: A Bayesian course with examples in R and Stan\" Secound Ed. 193-221.\n",
    "3. Vehtari Aki, Gelman Andrew, Gabry Jonah (2016) \"Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC\"\n"
   ]
  }
 ],
 "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.7.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}