Skip to content

Instantly share code, notes, and snippets.

@josef-pkt

josef-pkt/try_all_subset_regression2.ipynb

Created November 26, 2019 18:05
Show Gist options
  • Save josef-pkt/330f21972556852e1b0fc04bd8b3a3fc to your computer and use it in GitHub Desktop.
Save josef-pkt/330f21972556852e1b0fc04bd8b3a3fc to your computer and use it in GitHub Desktop.
Download ZIP
all subset regression and researcher degrees of freedom - sensitivity of parameter estimate to model specification
Display the source blob
Display the rendered blob
Raw
try_all_subset_regression2.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# All subsets regression with sweep algorithm"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from statsmodels.regression.linear_model import OLS\n",
"\n",
"from numpy.testing import assert_almost_equal\n",
"\n",
"from statsmodels.regression._stepwise import (\n",
" SequentialOLSQR, loglike_ssr, sweep, StepwiseOLSSweep, get_sweep_matrix_data,\n",
" all_subset)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"16384"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nobs, k_vars = 100, 15\n",
"\n",
"np.random.seed(85325783)\n",
"x = np.random.randn(nobs, k_vars)\n",
"x[:,0] = 1.\n",
"#x[:,1] = x[:,0] #to check what happens with singular exog matrix\n",
"beta = 1 / np.arange(1, k_vars + 1) **3\n",
"beta[3] += 0.5\n",
"y = x[:, :k_vars].dot(beta) + np.random.randn(nobs)\n",
"#last two variables have no effect, used to get aic, bic minimum in interior\n",
"\n",
"2**(k_vars - 1)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1.00000000e+00, 1.25000000e-01, 3.70370370e-02, 5.15625000e-01,\n",
" 8.00000000e-03, 4.62962963e-03, 2.91545190e-03, 1.95312500e-03,\n",
" 1.37174211e-03, 1.00000000e-03, 7.51314801e-04, 5.78703704e-04,\n",
" 4.55166136e-04, 3.64431487e-04, 2.96296296e-04])"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"beta"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ranking of explanatory variables by magnitude of beta parameter"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0, 3, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14],\n",
" dtype=int64)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.argsort(beta)[::-1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## All subset regression\n",
"\n",
"We use the `all_subset` function to compute the regression results for 8192 models. The running time for this is around 50 seconds.\n",
"\n",
"The `res_all` attribute contains the regression results in the sequence as the models where computed using the sweep algorithm.\n",
"The `sorted_frame` attribute contains a pandas DataFrame with results sorted by `aic`. From those results we can see that all models include the first two variables as required. Variable with index 3 is included in most modelsin the top 30 models ranked by aic.\n",
"The importance of explanatory variables in the simulated model depends on the beta parameter and explanatory variables sorted by importance are [0, 3, 1, 2, 4, 5]. Our sample size is small, only 100 observations, and the selection does not reflect the magnitude of the parameters in the model that was used to simulated the data with the exception of parameter with index 3."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.6455192565917969, 8192)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import time\n",
"t0 = time.time()\n",
"res_allf = all_subset(y, x, keep_exog=2)\n",
"t1 = time.time()\n",
"res_all0 = res_allf # keep for later\n",
"t1 - t0, len(res_allf.aic)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0 1] [152.8483953562869, 330.2153430543272]\n",
"[0 1 2] [152.84676549756713, 332.2142767249059]\n",
"[0 1 2 3] [107.31434995150929, 298.8469257849491]\n",
"[0 1 3] [107.71828611100258, 297.2226237645825]\n",
"[0 1 3 4] [107.66568239545632, 299.17377730961465]\n",
"[0 1 2 3 4] [107.25666440875483, 300.7931575308201]\n",
"[0 1 2 4] [152.84525430481443, 334.21328802212355]\n",
"[0 1 4] [152.84692459430477, 332.2143808138962]\n",
"[0 1 4 5] [152.77877164104606, 334.1697818451394]\n",
"[0 1 2 4 5] [152.7745384759872, 336.16701102582965]\n",
"[0 1 2 3 4 5] [106.85407272695764, 302.41709777628347]\n",
"[0 1 3 4 5] [107.19361093477245, 300.73435277927723]\n",
"[0 1 3 5] [107.26973789884704, 298.80534576496547]\n",
"[0 1 2 3 5] [106.93379840159686, 300.49168168781614]\n",
"[0 1 2 5] [152.77500037206934, 334.16731336376404]\n",
"[0 1 5] [152.7792146103344, 332.1700717863671]\n",
"[0 1 5 6] [147.15150059740927, 330.41695527885236]\n",
"[0 1 2 5 6] [146.93112591112796, 332.26708261881265]\n",
"[0 1 2 3 5 6] [105.69641733474626, 301.32778780659254]\n",
"[0 1 3 5 6] [105.82586699881458, 299.45018595545633]\n",
"[0 1 3 4 5 6] [105.74242875617813, 301.3713100098289]\n",
"[0 1 2 3 4 5 6] [105.6110419531209, 303.2469810176995]\n",
"[0 1 2 4 5 6] [146.9306991526969, 334.2667921704491]\n",
"[0 1 4 5 6] [147.15099641850796, 332.4166126525366]\n",
"[0 1 4 6] [147.22194566195608, 330.464816298103]\n",
"[0 1 2 4 6] [146.98868726928464, 332.30625068868164]\n",
"[0 1 2 3 4 6] [105.74981004466824, 301.3782902078115]\n",
"[0 1 3 4 6] [105.89264744886007, 299.51327014290626]\n",
"[0 1 3 6] [105.96262806891822, 297.57933469799684]\n",
"[0 1 2 3 6] [105.82217169245556, 299.44669402008833]\n"
]
}
],
"source": [
"#res_allf = all_subset(y, x, keep_exog=1)\n",
"for res_i in res_allf.res_all[:30]:\n",
" print(np.nonzero(res_i[1])[0], res_i[2:])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" exog_idx aic\n",
"0 [ 0 1 3 8 12] 290.938877\n",
"1 [ 0 1 3 6 8 12] 291.491601\n",
"2 [ 0 1 3 7 8 12] 292.390845\n",
"3 [ 0 1 2 3 8 12] 292.421399\n",
"4 [ 0 1 3 5 8 12] 292.491442\n",
"5 [ 0 1 3 8 11 12] 292.617453\n",
"6 [ 0 1 3 8 12 13] 292.672760\n",
"7 [ 0 1 3 8 10 12] 292.759138\n",
"8 [ 0 1 3 8 12 14] 292.897470\n",
"9 [ 0 1 3 8 9 12] 292.907810\n",
"10 [ 0 1 3 4 8 12] 292.938553\n",
"11 [ 0 1 3 6 8 11 12] 293.012504\n",
"12 [ 0 1 3 6 7 8 12] 293.025115\n",
"13 [ 0 1 3 6 8 12 14] 293.165864\n",
"14 [ 0 1 3 6 8 12 13] 293.232830\n",
"15 [ 0 1 2 3 6 8 12] 293.256914\n",
"16 [ 0 1 3 6 8 10 12] 293.278316\n",
"17 [ 0 1 3 5 6 8 12] 293.326357\n",
"18 [ 0 1 3 6 8 9 12] 293.445069\n",
"19 [ 0 1 3 4 6 8 12] 293.487363\n",
"20 [0 1 3 8] 293.801012\n",
"21 [ 0 1 2 3 7 8 12] 293.991639\n",
"22 [ 0 1 2 3 5 8 12] 294.049393\n",
"23 [ 0 1 2 3 8 12 13] 294.110370\n",
"24 [ 0 1 3 5 7 8 12] 294.118420\n",
"25 [ 0 1 3 7 8 11 12] 294.178963\n",
"26 [ 0 1 3 5 8 11 12] 294.228531\n",
"27 [ 0 1 2 3 8 11 12] 294.230661\n",
"28 [ 0 1 3 8 11 12 13] 294.237764\n",
"29 [ 0 1 3 7 8 12 13] 294.239461\n"
]
}
],
"source": [
"df_summ = res_allf.sorted_frame()\n",
"df_best = df_summ.iloc[:30]\n",
"print(df_best.to_string(formatters={'exog_idx': lambda x: str(x).ljust(30)}))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"exog_idx object\n",
"aic float64\n",
"dtype: object"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_summ.dtypes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## All subset regression with at most k_max variables\n",
"\n",
"In the following example we impose a restriction that we want to have at most k_max = 5 explanatory variables in a model. \n",
"\n",
"We can compute the number of candidate models using combinations. There are 4944 possible models, 3003 have 5 variables, if we don't require that some variables are contained in every model.\n",
"In this example, we want to always include the first two variables, i.e. the constant and one variable of interest. In this case we have 13 extra variables, and we want to choose at most 3 of them. So, we have 378 candidate models in this case."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(4944.0, [1.0, 15.0, 105.0, 455.0, 1365.0, 3003.0])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy.special import comb\n",
"k_ = 15\n",
"kmax_ = 5\n",
"comb(k_, np.arange(kmax_+1)).sum(), comb(k_, np.arange(kmax_+1)).tolist()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(378.0, [1.0, 13.0, 78.0, 286.0])"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k_ = 13\n",
"kmax_ = 3\n",
"comb(k_, np.arange(kmax_+1)).sum(), comb(k_, np.arange(kmax_+1)).tolist()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.29019641876220703, 377)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import time\n",
"t0 = time.time()\n",
"res_allf = all_subset(y, x, keep_exog=2, k_max=5)\n",
"t1 = time.time()\n",
"t1 - t0, len(res_allf.aic)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The all subset regression with the k_max restriction returns results for 377 models in 30 seconds.\n",
"\n",
"**BUG?** it looks like we are missing one model. This might be the model without any extra variables. (This is the empty model in terms of variables to select, and is currently not included in the candidate list.)\n",
"\n",
"We can again look at the sorted dataframe to see what our best models are. The best model from the unrestricted all subset regression has 5 explanatory variables, therefore it is also the best here. However, the second best model that has 5 or fewer explanatory variables was only the 21st best model in the unrestricted case. In this example, most of the top ranked models in terms of aic have 6 or 7 explanatory variables and do not satisfy our restriction on the maximum number of explanatory variables in a model."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" exog_idx aic\n",
"0 [ 0 1 3 8 12] 290.938877\n",
"1 [0 1 3 8] 293.801012\n",
"2 [0 1 3 6 8] 294.791141\n",
"3 [ 0 1 3 6 12] 295.270627\n",
"4 [ 0 1 3 8 13] 295.285584\n",
"5 [ 0 1 3 8 11] 295.334311\n",
"6 [0 1 2 3 8] 295.342805\n",
"7 [0 1 3 7 8] 295.381120\n",
"8 [0 1 3 8 9] 295.398978\n",
"9 [0 1 3 5 8] 295.464322\n",
"10 [ 0 1 3 12] 295.467308\n",
"11 [ 0 1 3 8 10] 295.567308\n",
"12 [ 0 1 3 8 14] 295.595339\n",
"13 [0 1 3 4 8] 295.800614\n",
"14 [ 0 1 3 7 12] 295.851438\n",
"15 [ 0 1 3 12 13] 296.753861\n",
"16 [ 0 1 3 11 12] 296.910295\n",
"17 [ 0 1 3 5 12] 296.935628\n",
"18 [ 0 1 2 3 12] 297.052866\n",
"19 [ 0 1 3 10 12] 297.219601\n"
]
}
],
"source": [
"df_summ = res_allf.sorted_frame()\n",
"df_best = df_summ.iloc[:20]\n",
"print(df_best.to_string(formatters={'exog_idx': lambda x: str(x).ljust(30)}))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can further restrict the number of explanatory variables so that only one extra variable is included. Here, the model with variable 3 has much lower aic than other models. Only models that include either variable three or six have lower aic than the \"empty\" model without extra regressors. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.08907365798950195 13\n",
" exog_idx aic\n",
"0 [0 1 3] 297.222624\n",
"1 [0 1 6] 328.465828\n",
"2 [0 1] 330.215343\n",
"3 [0 1 8] 330.346281\n",
"4 [ 0 1 12] 331.010611\n",
"5 [ 0 1 14] 331.030102\n",
"6 [ 0 1 13] 331.612151\n",
"7 [0 1 7] 331.633902\n",
"8 [0 1 9] 331.967544\n",
"9 [0 1 5] 332.170072\n",
"10 [ 0 1 11] 332.184214\n",
"11 [0 1 4] 332.214381\n",
"12 [ 0 1 10] 332.215240\n"
]
}
],
"source": [
"import time\n",
"t0 = time.time()\n",
"res_allf = all_subset(y, x, keep_exog=2, k_max=3)\n",
"t1 = time.time()\n",
"print(t1 - t0, len(res_allf.aic))\n",
"\n",
"df_summ = res_allf.sorted_frame()\n",
"df_best = df_summ.iloc[:20]\n",
"print(df_best.to_string(formatters={'exog_idx': lambda x: str(x).ljust(30)}))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we increase the number of extra variables to at most 2, then including either 8 or 12 in addition to 3 leads to a better aic than the model that only includes 3."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.14910459518432617 91\n",
" exog_idx aic\n",
"0 [0 1 3 8] 293.801012\n",
"1 [ 0 1 3 12] 295.467308\n",
"2 [0 1 3] 297.222624\n",
"3 [0 1 3 6] 297.579335\n",
"4 [0 1 3 7] 297.912861\n",
"5 [ 0 1 3 13] 298.234392\n",
"6 [ 0 1 3 11] 298.526839\n",
"7 [0 1 3 5] 298.805346\n",
"8 [0 1 2 3] 298.846926\n",
"9 [ 0 1 3 14] 298.867615\n",
"10 [ 0 1 3 10] 298.926273\n",
"11 [0 1 3 9] 299.023253\n",
"12 [0 1 3 4] 299.173777\n",
"13 [ 0 1 6 14] 327.676746\n",
"14 [0 1 6] 328.465828\n",
"15 [ 0 1 6 12] 328.720660\n",
"16 [0 1 6 8] 329.094295\n",
"17 [ 0 1 6 13] 329.891581\n",
"18 [0 1 6 7] 330.071934\n",
"19 [0 1 6 9] 330.199685\n"
]
}
],
"source": [
"import time\n",
"t0 = time.time()\n",
"res_allf = all_subset(y, x, keep_exog=2, k_max=4)\n",
"t1 = time.time()\n",
"print(t1 - t0, len(res_allf.aic))\n",
"\n",
"df_summ = res_allf.sorted_frame()\n",
"df_best = df_summ.iloc[:20]\n",
"print(df_best.to_string(formatters={'exog_idx': lambda x: str(x).ljust(30)}))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## TODO\n",
"\n",
"Currently there is no option to save additional results. We can store the parameters of interest and their standard deviation, `bse` in statsmodels, or their confidence intervals so that we can look at the sensitivity of the parameters of interest to this specification uncertainty.\n",
"\n",
"This example has uncorrelated explanatory variables. The more difficult case will be when there is high multicollinearity among the regressors."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sensitivity of parameter of interest to model specification uncertainty\n",
"\n",
"I changed the code so that the parameters of interest, those that are required to be included in the model. are also stored for each subset regression.\n",
"\n",
"Below are some experiments for how to look at the results for the parameter 1, i.e. the first slope parameter."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"params = np.asarray(res_all0.params_keep_all)\n",
"bse = np.asarray(res_all0.bse_keep_all)\n",
"df_resid = np.asarray(res_all0.df_resid_all)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[array([0.9846275 , 0.06720861]),\n",
" array([0.98474711, 0.06732363]),\n",
" array([1.02187845, 0.08201408]),\n",
" array([1.02340181, 0.08367705]),\n",
" array([1.02469827, 0.08590344]),\n",
" array([1.02322651, 0.08433521]),\n",
" array([0.98453795, 0.06695176]),\n",
" array([0.98441971, 0.06684038]),\n",
" array([0.98572556, 0.06508744]),\n",
" array([0.98594454, 0.06522531])]"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res_all0.params_keep_all[:10]"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[array([0.1248947 , 0.12610032]),\n",
" array([0.12559125, 0.12679842]),\n",
" array([0.10594154, 0.10682323]),\n",
" array([0.10556199, 0.10643528]),\n",
" array([0.10625322, 0.10745494]),\n",
" array([0.10663594, 0.10784511]),\n",
" array([0.12642548, 0.12802675]),\n",
" array([0.12572035, 0.12731985]),\n",
" array([0.12650275, 0.1282328 ]),\n",
" array([0.12723666, 0.12893216])]"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res_all0.bse_keep_all[:10]"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"idx_sort_params1 = np.argsort(params[:,1])\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0xf828dca508>]"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAZOUlEQVR4nO3dfbBU9X3H8ffHS8GH+IAGnQS4hUTahIyOmhVMba0TS0TT5voHKVdNim0aJu04k8S0HTI4TaWx0bZjSadMG3zIxIeIkTbkTiRhqMZOJ1HKEhVEJV4JgQs2kqAmYx4I8O0fe65dl73cs3fP7tk9+3nNMPc8/M7Z3zlz+Nzf/Z2z56eIwMzMiuu4vCtgZmat5aA3Mys4B72ZWcE56M3MCs5Bb2ZWcJPyrkCtN7/5zTFr1qy8q2Fm1lW2bNnyo4iYVm9dxwX9rFmzKJfLeVfDzKyrSPrBWOvcdWNmVnAOejOzgnPQm5kVnIPezKzgHPRmZgXXcU/dmJn1mhvXbeP+TXs4HEGfxNXzZ/LZq87JbP8OejOzHN24bhv3Pr779fnDEa/PZxX27roxM8vR/Zv2NLR8Ihz0ZmY5OjzGmCBjLZ8IB72ZWcE56M3MCs5Bb2ZWcA56M7OCc9CbmRWcg97MLCfvWL6+LZ/joDczy8kvDmf3COWxOOjNzAouVdBLWihph6RhScvqrL9E0nclHZK0qGr5r0vaIulJSdslfSzLypuZ2fjGfdeNpD5gFbAAGAE2SxqKiGeqiu0GrgP+ombzF4HfiohfSnoT8HSy7b5Mam9mVlArF5+X2b7SvNRsHjAcETsBJK0BBoDXgz4idiXrjlRvGBEHq2an4K4iMzNg/BuxV50/PbPPShO804Hqt+uMJMtSkTRT0tZkH7fWa81LWiqpLKm8f//+tLs2M+ta7boRC+mCXnWWpa5hROyJiHOBs4Elks6qU2Z1RJQiojRt2rS0uzYz60rrntjb1s9LE/QjwMyq+RlAw33sSUt+O/A7jW5rZlYkn3jgybZ+Xpqg3wzMkTRb0mRgEBhKs3NJMySdkExPBS4Gdky0smZm3W7+zRvHLZPljVhIEfQRcQi4HtgAPAt8JSK2S1oh6QMAki6UNAJ8EPiCpO3J5u8ENkl6Cvgv4B8jYlumR2Bm1kV++NOD45bJ8kYspBxKMCLWA+trlv111fRmKl06tdttBM5tso5mZoVw7me+OW6ZSfXuijbJjzuambXBtbc/xk9+eXjccsOfe3/mn+3Bwc3MWmzWsody/Xy36M3MWqiRkN91S/ateXDQm5m1TCMhf3xfCzrnEw56M7MWaLS75rmbr2xRTdxHb2aWqYn0x7eqy2aUg97MLAM3rtvGvY/vbni7Voc8OOjNzJqy4LZHef6l1ya0bTtCHhz0ZmYT8o7l65t6A2W7Qh4c9GZmDWv2ufh2hjw46M3MUsnqS0/tDnlw0JuZHdNEb7LWWrn4vMxfVpaWg97MrI5m++Cr5dGKr+agNzOrkuV7afIO+FEOejMzihnwo1IFvaSFwOeBPuCOiLilZv0lwEoq754fjIi1yfLzgH8FTgEOAzdHxAPZVd/MbOKyfqtkpwX8qHGDXlIfsApYQGX82M2ShiLimapiu4HrgL+o2fxnwB9FxPOS3gpskbQhIl7JpPZmZg1q5gtOY+nUgB+VpkU/DxiOiJ0AktYAA8DrQR8Ru5J1R6o3jIjvVU3vk/QSMA1w0JtZ26x7Ym9LBuTu9IAflSbopwN7quZHgPmNfpCkecBk4IU665YCSwH6+/sb3bWZWV2tGvCjWwJ+VJqgr/eS5IaeOZL0FuAeYElEHKldHxGrgdUApVIpm+eZzKwntXI0p24L+FFpgn4EmFk1PwPYl/YDJJ0CPATcGBGPN1Y9M7PxtTLc55x5EhtvuLRl+2+HNEG/GZgjaTawFxgErkmzc0mTga8Cd0fEgxOupZlZjXM/881Ug21P1MVvP537Pvqelu2/ncYN+og4JOl6YAOVxyvviojtklYA5YgYknQhlUCfCvyBpJsi4l3AHwKXAGdIui7Z5XURkf1dETPrCa0eaLtbu2eORRGd1SVeKpWiXC7nXQ0z6yAO9/FJ2hIRpXrr/M1YM+tYvrGaDQe9mXWcVgV8nm+QzJOD3sw6QpZvi6zVS633ehz0ZpYrf6mp9Rz0ZpaLrAb0qOZwr89Bb2ZtV+RXAnciB72ZtVU3j73arRz0ZtY2zYb8hy7q57NXnZNRbXqHg97MWq7Z1wS79d4cB72ZtVQzA3044LPhoDezlrn29scmFPIO+Gw56M2sZb79woGGyjvgW+O4vCtgZsXU6I1Xh3zrOOjNLFfH98kh32LuujGzzKVtzTvg28MtejPLhUO+fVIFvaSFknZIGpa0rM76SyR9V9IhSYtq1n1T0iuSvp5Vpc2su1389tPzrkJPGTfoJfUBq4ArgLnA1ZLm1hTbDVwHfLnOLv4B+HBz1TSzIinKWKzdIk2Lfh4wHBE7I+IgsAYYqC4QEbsiYitwpHbjiHgY+GkWlTWzzrfgtkePud5dNu2XJuinA3uq5keSZZmRtFRSWVJ5//79We7azNpsot+CtdZJE/SqsyzTYWAiYnVElCKiNG3atCx3bWbW89IE/Qgws2p+BrCvNdUxM7OspQn6zcAcSbMlTQYGgaHWVsvMishP2+Rj3KCPiEPA9cAG4FngKxGxXdIKSR8AkHShpBHgg8AXJG0f3V7SfwMPApdJGpF0eSsOxMw6n5+2yUeqb8ZGxHpgfc2yv66a3kylS6fetr/TTAXNzKw5/masmWXqQxf1N7TcWs/vujGzTI0O9Xf/pj0cjqBP4ur5Mz0EYI4UkemTkk0rlUpRLpfzroaZWVeRtCUiSvXWuevGzKzgHPRmZgXnoDczKzjfjDWzzK17Yi//sGEH+175OW897QT+8vLf5KrzM31FljXAQW9mmVr3xF4+9eBTHD5SedBj7ys/51MPPgXgsM+Ju27MLFPLv7rt9ZAfdfhIsPyr23KqkTnozSxTrx083NByaz0HvZlZwTnozcwKzkFvZlZwDnozy8x448VaPhz0ZpYZjxfbmVIFvaSFknZIGpa0rM76SyR9V9IhSYtq1i2R9Hzyb0lWFTczs3TGDXpJfcAq4ApgLnC1pLk1xXYD1wFfrtn2dOAzwHxgHvAZSVObr7aZdZpZyx7Kuwo2hjQt+nnAcETsjIiDwBpgoLpAROyKiK3AkZptLwc2RsSBiHgZ2AgszKDeZtZB0oT8ysXntaEmVk+aoJ8O7KmaH0mWpdHMtmbWBdK25P36g/ykedeN6ixLO1pJqm0lLQWWAvT3e7gxs25w47pt3Pv47ryrYSmkCfoRYGbV/AxgX8r9jwCX1mz7aG2hiFgNrIbKCFMp921mOWm0P/6UKX0tqomlkSboNwNzJM0G9gKDwDUp978B+LuqG7DvAz7dcC3NrCNM9Ibr1pt8ay5P4wZ9RBySdD2V0O4D7oqI7ZJWAOWIGJJ0IfBVYCrwB5Juioh3RcQBSX9L5ZcFwIqIONCiYzGzFmnmiZpdt7w/w5rYRHhwcDOra90Te/nEA082tQ+HfPsca3BwDzxiZm8we9lDqZ+2OBaHfOdw0JsZkO0XnhzyncVBb9bDsv42qwO+MznozXrQgtsezfQFZMf3ieduvjKz/Vm2HPRmPSTrFvxZJ09m0/IFme7TsuegN+sBWQf8ysXn+ZUGXcRBb1ZgWXfRuA++OznozQroHcvX84vD2X1HxgHf3Rz0ZgXip2isHge9WUFkFfIO9+Jx0Jt1sWtvf4xvv5DN66Mc8MXloDfrUm7BW1oOerMukuVgHw743uGgN+sCfg+NNcNBb9bBsnxM0gHfuxz0Zh3IXTSWpVRBL2kh8HkqI0zdERG31KyfAtwNvBv4MbA4InZJmgx8ASgBR4CPR8Sj2VXfrHh8k9WyNm7QS+oDVgELqAz2vVnSUEQ8U1XsI8DLEXG2pEHgVmAx8FGAiDhH0pnANyRdGBFHsj4Qs26XVSveAW+10rTo5wHDEbETQNIaYACoDvoB4G+S6bXAv0gSMBd4GCAiXpL0CpXW/f9kUnuzgsiiFe8XjdlY0gT9dGBP1fwIMH+sMslg4q8CZwBPAQPJL4eZVLp2ZlIT9JKWAksB+vv7Gz8Ksy6U1QvH3IK38aQJetVZVvsYwFhl7gLeCZSBHwDfAQ4dVTBiNbAaKoODp6iTWddyF421W5qgH6HSCh81A9g3RpkRSZOAU4EDERHAJ0cLSfoO8HxTNTbrUr7JanlJE/SbgTmSZgN7gUHgmpoyQ8AS4DFgEfBIRISkEwFFxGuSFgCHam7imhWW3yRpnWLcoE/63K8HNlB5vPKuiNguaQVQjogh4E7gHknDwAEqvwwAzgQ2SDpC5ZfEh1txEGadIutwh0q/6Pcd8tYEVXpXOkepVIpyuZx3NcxSmX/zRn7404Mt279b8ZaWpC0RUaq3zt+MNWtAK1rs9fhRScuSg97sGNoV7KPcgrdWcNCbVWl3sI9ywFsrOeit5+UV7r7Jau3ioLeek1ewg1vulg8HvRVanqEObrVbZ3DQW2HkHeqj3Gq3TuOgt67UKaF+8dtP576Pvifvapgdk4PeukKWQ+o1w8+3Wzdy0FtH6pQWO7grxrqfg946RqeEu4PdisZBb7nplGD/0EX9fPaqc/KuhlnLOOitrfIOdz/uaL3IQW8tde3tj/HtFw7k8tluqZtVOOgtc35fjFlncdBbJvIIdwe7WTqpgl7SQuDzVEaYuiMibqlZPwW4G3g38GNgcUTskvRrwB3ABcln3R0Rn8uw/pYTB7tZ9xg36CX1AauABVQGAd8saahm7NePAC9HxNmSBoFbgcXAB4EpEXFOMn7sM5Luj4hdWR+Itd66J/byiQeebNvnzTnzJDbecGnbPs+sqNK06OcBwxGxE0DSGmAAqA76AeBvkum1wL9IEhDASZImAScAB4GfZFN1a4cFtz3K8y+91pbPmiQY/pxb7WZZSxP004E9VfMjwPyxyiSDib8KnEEl9AeAF4ETgU9GxFGPYEhaCiwF6O/vb/AQrBXa1TXjd8WYtV6aoFedZbUvHRmrzDzgMPBWYCrw35L+c/Svg9cLRqwGVkNlcPAUdbIWaUfAu6/drL3SBP0IMLNqfgawb4wyI0k3zanAAeAa4JsR8SvgJUnfBkrATqyjtDrgHe5m+UkT9JuBOZJmA3uBQSoBXm0IWAI8BiwCHomIkLQbeK+ke6l03VwErMyq8ta8G9dt497Hd2e+X7/l0axzjBv0SZ/79cAGKo9X3hUR2yWtAMoRMQTcCdwjaZhKS34w2XwV8EXgaSrdO1+MiK0tOA6bgCxb8b6Rata5FNFZXeKlUinK5XLe1Si8rELeXTJmnUHSlogo1Vvnb8b2oGZD3uFu1l0c9D2kmVGazjp5MpuWL8i4RmbWDg76HjHRVryfczfrfg76HjCRkPfNVbPicNAX3LW3P9bwNu6DNysWB32BNfqeGo++ZFZMDvqCarS7xl9wMisuB33BTOSbru6qMSs2B32BTOSmq0PerPgc9AUw0UcnHfJmveG4vCtgzZloyK9cfF7GNTGzTuUWfZdq5jUGvvFq1lsc9F2m2XFb3V1j1nsc9F3ELyMzs4lw0HcBB7yZNSPVzVhJCyXtkDQsaVmd9VMkPZCs3yRpVrL8WklPVv07Isl3AVOateyhpkL+4ref7pA3s/Fb9JL6qIwUtYDK2LCbJQ1FxDNVxT4CvBwRZ0saBG4FFkfEfcB9yX7OAb4WERPvYO4Rjb66oB4HvJmNStOinwcMR8TOiDgIrAEGasoMAF9KptcCl0lSTZmrgfubqWwvmLXsoaZC/qyTJzvkzewN0vTRTwf2VM2PAPPHKpOMMfsqcAbwo6oyizn6F4QlshjazwFvZvWkCfraljlA7TBFxywjaT7ws4h4uu4HSEuBpQD9/f0pqlQcDngza7U0QT8CzKyanwHsG6PMiKRJwKnAgar1gxyj2yYiVgOroTI4eIo6dT0Pzm1m7ZIm6DcDcyTNBvZSCe1rasoMAUuAx4BFwCMREQCSjgM+CFySVaW7WRY3Wj36k5k1YtygT/rcrwc2AH3AXRGxXdIKoBwRQ8CdwD2Shqm05AerdnEJMBIRO7Ovfvdo9huto9yCN7NGKWl4d4xSqRTlcjnvamTGXTRm1g6StkREqd46fzO2RbIK+FOm9LH1poWZ7MvMepODvgWyCPk5Z57Exhsubb4yZtbzHPQZyiLgPUC3mWXNQZ8B98ObWSdz0DfBAW9m3cBB36Cswh0c8GbWHg76lBzwZtatHPTHkGW4gwPezPLhoK+RdbiPcsibWV4c9Iks3kFTjwPezPLmoKc1rXgHvJl1ip4P+ixD/vg+8dzNV2a2PzOzLPR00Ps5eDPrBT0Z9A54M+slPRX0HrbPzHpRzwR9MyHvcDezbpYq6CUtBD5PZYSpOyLilpr1U4C7gXcDPwYWR8SuZN25wBeAU4AjwIUR8YusDiCNiYa8A97MiuC48QpI6gNWAVcAc4GrJc2tKfYR4OWIOBv4J+DWZNtJwL3AxyLiXcClwK8yq30KDnkz63XjBj0wDxiOiJ0RcRBYAwzUlBkAvpRMrwUukyTgfcDWiHgKICJ+HBGHs6n6+BzyZmbpum6mA3uq5keA+WOVSQYTfxU4A/gNICRtAKYBayLi72s/QNJSYClAf39/o8dQ10RC3gFvZkWUJuhVZ1ntiOJjlZkE/DZwIfAz4OFkANuH31AwYjWwGiqDg6eo0zHNbjDkVy4+j6vOn97sx5qZdaQ0QT8CzKyanwHsG6PMSNIvfypwIFn+XxHxIwBJ64ELgIdpkXcsX3/Ub6FjcSvezIouTR/9ZmCOpNmSJgODwFBNmSFgSTK9CHgkIgLYAJwr6cTkF8DvAs9kU/WjLbjtUX5xOF3Mn3XyZIe8mfWEcVv0SZ/79VRCuw+4KyK2S1oBlCNiCLgTuEfSMJWW/GCy7cuSbqPyyyKA9RHRmvcAQ0Nvn9y0fEGrqmFm1lFUaXh3jlKpFOVyueHtGrn56pa8mRVNcv+zVG9dmq6bQnHIm1mv6amgd8ibWS8qRNAvuO3Rccs45M2sVxUi6Me7CTvnzJPaVBMzs85TiKAfz8YbLs27CmZmuemJoDcz62WFD/qzTp6cdxXMzHJV+KD3F6PMrNcVPujNzHqdg97MrOAKEfSnTOlraLmZWS8pRNBvvWnhUaF+ypQ+tt60MKcamZl1jlSDg3cDh7qZWX2FaNGbmdnYHPRmZgXnoDczKzgHvZlZwTnozcwKruOGEpS0H/hBE7t4M/CjjKpTZD5P6fg8pedzlU6rztOvR8S0eis6LuibJak81riJ9v98ntLxeUrP5yqdPM6Tu27MzArOQW9mVnBFDPrVeVegS/g8pePzlJ7PVTptP0+F66M3M7M3KmKL3szMqjjozcwKrjBBL2mhpB2ShiUty7s+7SZppqRvSXpW0nZJH0+Wny5po6Tnk59Tk+WS9M/J+doq6YKqfS1Jyj8vaUlex9RKkvokPSHp68n8bEmbkmN+QNLkZPmUZH44WT+rah+fTpbvkHR5PkfSWpJOk7RW0nPJtfUeX1NHk/TJ5P/d05Lul3R8R11TEdH1/4A+4AXgbcBk4Clgbt71avM5eAtwQTJ9MvA9YC7w98CyZPky4NZk+krgG4CAi4BNyfLTgZ3Jz6nJ9NS8j68F5+sG4MvA15P5rwCDyfS/AX+WTP858G/J9CDwQDI9N7nOpgCzk+uvL+/jasF5+hLwp8n0ZOA0X1NHnaPpwPeBE6qupes66ZoqSot+HjAcETsj4iCwBhjIuU5tFREvRsR3k+mfAs9SuQAHqPxnJfl5VTI9ANwdFY8Dp0l6C3A5sDEiDkTEy8BGoFAv+5c0A3g/cEcyL+C9wNqkSO15Gj1/a4HLkvIDwJqI+GVEfB8YpnIdFoakU4BLgDsBIuJgRLyCr6l6JgEnSJoEnAi8SAddU0UJ+unAnqr5kWRZT0r+FDwf2AScFREvQuWXAXBmUmysc9YL53Il8FfAkWT+DOCViDiUzFcf8+vnI1n/alK+F87T24D9wBeTbq47JJ2Er6k3iIi9wD8Cu6kE/KvAFjromipK0KvOsp58blTSm4B/Bz4RET85VtE6y+IYywtB0u8DL0XElurFdYrGOOsKfZ4Sk4ALgH+NiPOB16h01YylJ89Vco9igEp3y1uBk4Ar6hTN7ZoqStCPADOr5mcA+3KqS24k/RqVkL8vIv4jWfzD5M9nkp8vJcvHOmdFP5cXAx+QtItKF997qbTwT0v+7IY3HvPr5yNZfypwgOKfJ6gc40hEbErm11IJfl9Tb/R7wPcjYn9E/Ar4D+C36KBrqihBvxmYk9zlnkzlBsdQznVqq6SP707g2Yi4rWrVEDD6lMMS4GtVy/8oeVLiIuDV5M/wDcD7JE1NWirvS5YVQkR8OiJmRMQsKtfJIxFxLfAtYFFSrPY8jZ6/RUn5SJYPJk9QzAbmAP/TpsNoi4j4X2CPpN9MFl0GPIOvqVq7gYsknZj8Pxw9T51zTeV9xzrDO99XUnnS5AVged71yeH4f5vKn3lbgSeTf1dS6ft7GHg++Xl6Ul7AquR8bQNKVfv6Eyo3goaBP8772Fp4zi7l/5+6eVvyn2oYeBCYkiw/PpkfTta/rWr75cn52wFckffxtOgcnQeUk+tqHZWnZnxNHX2ebgKeA54G7qHy5EzHXFN+BYKZWcEVpevGzMzG4KA3Mys4B72ZWcE56M3MCs5Bb2ZWcA56M7OCc9CbmRXc/wEbUH15HumE/wAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(params[idx_sort_params1, 1], 'o')"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0xf829306888>]"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAATcUlEQVR4nO3df6zddX3H8efbS4sXjBahLvYW1hqbum5s4q7VzYU4GLa4zRKDs7goJGT9Z2TOuC4lW2JkWdR0m7iEGInCEKegjLHGJTYGTLZgYL21SilQqQzh3rJRAxeNNukP3vvjfA87HM7p/d7ec+4593Oej6Tp+f4453w++bav+7nv7+f7/UZmIkkq16sG3QBJUn8Z9JJUOINekgpn0EtS4Qx6SSrcGYNuQLvzzjsv16xZM+hmSNKSsnfv3p9k5spO24Yu6NesWcPU1NSgmyFJS0pE/LjbNks3klQ4g16SCmfQS1LhDHpJKpxBL0mFG7pZN5I0au7ZN8PO3Qc5PHuUVSvG2b5pPVdcNNGzzzfoJWmA7tk3w/V37+fo8ZMAzMwe5fq79wP0LOwt3UjSAO3cffClkG86evwkO3cf7Nl3OKKXpAFolmtmZo923H64y/rTYdBL0iJrL9d0smrFeM++z9KNJC2yTuWaVuPLxti+aX3Pvs8RvSQtkrnKNQATzrqRpKWlNdwDONVTuidWjHP/jkt63gaDXpL6pL0Wf6qQ73W5ppVBL0l9Mlctvqkf5ZpWBr0k9VidWnxTv8o1rQx6SeqhOlMnm/pZrmll0EtSD81VrmmekO13uaaVQS9JPXSqK1oXM9xbGfSSdJpa7zr5uvFlRHSfWbMYtfhuDHpJOg3ttfjZo8e77rtYtfhuDHpJOg3DMnWyDoNektp0KsnM/uI4q1aM87tvWcl3HjtSa+pkwMDKNa0Mekmi+60KWksyM7NH+coDT9X+zF7egXIhDHpJI28+tyqoa9B1+VYGvaSR0qks8/wvup9IrWtFW4ln0HX5Vga9pCJ1C/RuZZmFGOTUyToMeknFqFNn70VZptUwlWi6MeglLQndZsJ0G60vNNCbn9VekmnOujk8e3ToSjTdGPSSht6pLk7qx2h9GOa+95JBL2mo9OtkaR3jy8b41PsvLCbgmwx6SYtirtJL8/XPj53g+MnG2LxXJ0vbdSvLlDSKb2XQS+qp+c526fa6F0Yt0Lsx6CUt2CBmu3QziPu9D7taQR8Rm4HPAWPAFzPz023bLwZuBH4d2JqZd1Xr3wp8HngtcBL428y8s3fNl7SY6ozWFyPQV3Qp/YziaL2OOYM+IsaAm4DLgGlgT0TsysxHWnZ7CrgG+Iu2t/8C+EhmPh4Rq4C9EbE7M2d70npJfTHf8stijdZLPVnab3VG9BuBQ5n5BEBE3AFsAV4K+sx8str2YusbM/OHLa8PR8SzwErAoJcGZL7z0Rc70Je9KnjNq89wtN5DdYJ+Ani6ZXkaeMd8vygiNgLLgR912LYN2AZwwQUXzPejJc2hTg19MQK908lRw7z/6gR9dFg3r38HEfFG4Hbg6sx8sX17Zt4M3AwwOTm5WL8FSkOt7nTEOq9bpywuxn8wZ7sMlzpBPw2c37K8Gjhc9wsi4rXAvwN/nZkPzK950miZ78h7vq/7ydkuw6tO0O8B1kXEWmAG2Ap8qM6HR8Ry4F+BL2fmN067ldKQ68XoexCzV06Ho/WlZ86gz8wTEXEdsJvG9MpbMvNARNwATGXmroh4O41APwf4w4j4ZGb+KvBHwMXAuRFxTfWR12Tm9/vRGanXFno153xfD1O4G+jliMxh+qfVqNFPTU0NuhkaAfOdfVIST4qWJyL2ZuZkp21eGaviLeSS/FJCvnXKoiE+egx6LQmnWwMf9JzwhZpr5F3ntcEug15D41RhvpAa+FII9FbOXlGvGfTqm/mMwuuWUoZVt6s5HX1rGBj0Oi0LuYx+qY28PXGppc6g1yv0KsSHNbhbGeIaBQa9gOG5F0q/OCdco8ygHyGtI/XWp9m3h/swhvhCauAGukadQV+4biP1mdmjfOWBp17ab1jC3ZG31HsG/RI231r6IMJ8PvPADXSpPwz6IdKP6Yj9DHdPZEpLg0E/YAu5La0hLqkOg36A7tk3w/V37+fo8ZPA4OvkXpEplcmgH6Cduw++FPKLpTXMm7NumrNwDHepTAb9ALSWaxaDI3VptBn0i6y9XLNQ1tIlzcWgXyR1R/FOR5TUawb9Iqg7ire0IqkfDPo+mk8tfmLFOPfvuGQRWiVp1Bj0PdR+wVPrwzJOZXzZGNs3rV+EFkoaRQb9AtW54OlULNdI6jeD/jR0C/f5XPA0vmyMT73/QgNeUt8Z9DX1ItybHMVLWkwGfQ29ulWBo3hJg2DQn8JCr2BtfViGc94lDYpB38XpXsHq7QYkDRuDvov53HDMcJc0zAz6Lg7XvFWB4S5p2Bn0XaxaMd61Nm+4S1pKXjXoBgyr7ZvWM75s7GXrxpeNceMH38r9Oy4x5CUtGY7ou2gGefOWBs6akbRUGfSncMVFEwa7pCXP0o0kFa5W0EfE5og4GBGHImJHh+0XR8T3IuJERFzZtu1bETEbEd/sVaMlSfXNGfQRMQbcBFwObACuiogNbbs9BVwDfLXDR+wEPrywZkqSTledEf1G4FBmPpGZx4A7gC2tO2Tmk5n5EPBi+5sz817gZ71orCRp/uoE/QTwdMvydLWuZyJiW0RMRcTUkSNHevnRkjTy6gR9dFh3ujdw7Cgzb87MycycXLlyZS8/WpJGXp2gnwbOb1leDRzuT3MkSb1WJ+j3AOsiYm1ELAe2Arv62yxJUq/MGfSZeQK4DtgNPAp8PTMPRMQNEfE+gIh4e0RMAx8AvhARB5rvj4j/BL4BXBoR0xGxqR8dkSR1Fpk9Lbcv2OTkZE5NTQ26GZK0pETE3syc7LTNK2MlqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwtYI+IjZHxMGIOBQROzpsvzgivhcRJyLiyrZtV0fE49Wfq3vVcElSPXMGfUSMATcBlwMbgKsiYkPbbk8B1wBfbXvv64FPAO8ANgKfiIhzFt5sSVJddUb0G4FDmflEZh4D7gC2tO6QmU9m5kPAi23v3QR8OzOfy8zngW8Dm3vQbklSTXWCfgJ4umV5ulpXR633RsS2iJiKiKkjR47U/GhJUh11gj46rMuan1/rvZl5c2ZOZubkypUra360JKmOOkE/DZzfsrwaOFzz8xfyXklSD9QJ+j3AuohYGxHLga3Arpqfvxt4T0ScU52EfU+1TpK0SOYM+sw8AVxHI6AfBb6emQci4oaIeB9ARLw9IqaBDwBfiIgD1XufA/6Gxg+LPcAN1TpJ0iKJzLrl9sUxOTmZU1NTg26GJC0pEbE3Myc7bfPKWEkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKlytoI+IzRFxMCIORcSODtvPjIg7q+0PRsSaav3yiLg1IvZHxA8i4t09bb0kaU5zBn1EjAE3AZcDG4CrImJD227XAs9n5puBzwKfqdb/CUBmXghcBvx9RPhbhCQtojqhuxE4lJlPZOYx4A5gS9s+W4Dbqtd3AZdGRND4wXAvQGY+C8wCk71ouCSpnjpBPwE83bI8Xa3ruE9mngBeAM4FfgBsiYgzImIt8JvA+e1fEBHbImIqIqaOHDky/15IkrqqE/TRYV3W3OcWGj8YpoAbge8CJ16xY+bNmTmZmZMrV66s0SRJUl1n1NhnmpePwlcDh7vsMx0RZwCvA57LzAQ+1twpIr4LPL6gFkuS5qXOiH4PsC4i1kbEcmArsKttn13A1dXrK4H7MjMj4qyIOBsgIi4DTmTmIz1quySphjlH9Jl5IiKuA3YDY8AtmXkgIm4ApjJzF/Al4PaIOAQ8R+OHAcAbgN0R8SIwA3y4H52QJHUXjerK8JicnMypqalBN0OSlpSI2JuZHWc1Oqddkgpn0EtS4Qx6SSpcnemVI+WefTPs3H2Qw7NHWbVinO2b1nPFRe3Xh0nS0mHQt7hn3wzX372fo8dPAjAze5Tr794PYNhLWrIs3bTYufvgSyHfdPT4SXbuPjigFknSwhn0LQ7PHp3XeklaCgz6FqtWjM9rvSQtBQZ9i+2b1jO+bOxl68aXjbF90/oBtUiSFs6TsS2aJ1yddSOpJAZ9mysumjDYJRXF0o0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOK2PxYSOSyjbyQe/DRiSVbmSDvjmKn+lwr/nmw0YMekklGMmgbx/Fd+LDRiSVYiRPxnZ6ZGA7HzYiqRQjNaI/VbmmlQ8bkVSSkQn6OuUagAln3UgqzMgE/VzlmvFlY3zq/Rca8JKKMzJBf6qTq47iJZWs+KBv1uWzy/aJFePcv+OSRW2TJC2mooN+rrq8J10ljYKig/5UdXnLNZJGRdFB360uH2C5RtLIqHXBVERsjoiDEXEoInZ02H5mRNxZbX8wItZU65dFxG0RsT8iHo2I63vb/M7u2TfDuz59X9e6vBdDSRolc47oI2IMuAm4DJgG9kTErsx8pGW3a4HnM/PNEbEV+AzwQeADwJmZeWFEnAU8EhFfy8wne92R1ouhArqGvHV5SaOmzoh+I3AoM5/IzGPAHcCWtn22ALdVr+8CLo2IZt6eHRFnAOPAMeCnPWl5i+ZJ1+YVr6eaYeNceUmjpk6NfgJ4umV5GnhHt30y80REvACcSyP0twDPAGcBH8vM59q/ICK2AdsALrjggnl2od69a6zLSxpVdUb00WFd+6C52z4bgZPAKmAt8PGIeNMrdsy8OTMnM3Ny5cqVNZr0cnXuNGldXtKoqhP008D5LcurgcPd9qnKNK8DngM+BHwrM49n5rPA/cDkQhvdbq4Qty4vaZTVCfo9wLqIWBsRy4GtwK62fXYBV1evrwTuy8wEngIuiYazgXcCj/Wm6f9v+6b1jC8be9m65q8Y1uUljbo5a/RVzf06YDcwBtySmQci4gZgKjN3AV8Cbo+IQzRG8lurt98E3Ao8TCN7b83Mh3rdiWaI+9xXSXqlaAy8h8fk5GROTU0NuhmStKRExN7M7FgaH8knTEnSKDHoJalwBr0kFc6gl6TCGfSSVLihm3UTEUeAHy/gI84DftKj5iwVo9hnGM1+j2KfYTT7Pd8+/3Jmdry1wNAF/UJFxFS3KUalGsU+w2j2exT7DKPZ71722dKNJBXOoJekwpUY9DcPugEDMIp9htHs9yj2GUaz3z3rc3E1eknSy5U4opcktTDoJalwxQR9RGyOiIMRcSgidgy6Pf0SEedHxHci4tGIOBARH63Wvz4ivh0Rj1d/nzPotvZaRIxFxL6I+Ga1vDYiHqz6fGf1vISiRMSKiLgrIh6rjvlvlX6sI+Jj1b/thyPiaxHx6hKPdUTcEhHPRsTDLes6HtvqmR7/WOXbQxHxtvl8VxFBHxFjNO59fzmwAbgqIjYMtlV9cwL4eGb+Co0Hufxp1dcdwL2ZuQ64t1ouzUeBR1uWPwN8turz88C1A2lVf32OxlPa3gL8Bo3+F3usI2IC+DNgMjN/jcYzMLZS5rH+J2Bz27pux/ZyYF31Zxvw+fl8URFBT+PZtIcy84nMPAbcQeOh5MXJzGcy83vV65/R+I8/QaO/t1W73QZcMZgW9kdErAZ+H/hitRzAJTQeQA9l9vm1wMU0HuxDZh7LzFkKP9Y0Hog0Xj2W9CzgGQo81pn5HzQe1NSq27HdAnw5Gx4AVkTEG+t+VylBPwE83bI8Xa0rWkSsAS4CHgR+KTOfgcYPA+ANg2tZX9wI/CXwYrV8LjCbmSeq5RKP+ZuAI8CtVcnqi9UjOYs91pk5A/wdjceQPgO8AOyl/GPd1O3YLijjSgn66LCu6HmjEfEa4F+AP8/Mnw66Pf0UEX8APJuZe1tXd9i1tGN+BvA24POZeRHwcwoq03RS1aS3AGuBVcDZNMoW7Uo71nNZ0L/3UoJ+Gji/ZXk1cHhAbem7iFhGI+T/OTPvrlb/b/NXuervZwfVvj54F/C+iHiSRlnuEhoj/BXVr/dQ5jGfBqYz88Fq+S4awV/ysf494L8z80hmHgfuBn6b8o91U7dju6CMKyXo9wDrqjPzy2mcvNk14Db1RVWb/hLwaGb+Q8umXcDV1eurgX9b7Lb1S2Zen5mrM3MNjWN7X2b+MfAd4Mpqt6L6DJCZ/wM8HRHrq1WXAo9Q8LGmUbJ5Z0ScVf1bb/a56GPdotux3QV8pJp9807ghWaJp5bMLOIP8F7gh8CPgL8adHv62M/fofEr20PA96s/76VRs74XeLz6+/WDbmuf+v9u4JvV6zcB/wUcAr4BnDno9vWhv28FpqrjfQ9wTunHGvgk8BjwMHA7cGaJxxr4Go3zEMdpjNiv7XZsaZRubqrybT+NWUm1v8tbIEhS4Uop3UiSujDoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuH+D/iy1A9hEYXyAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"idx_sort_aic = np.argsort(res_all0.aic)\n",
"k_best = 100\n",
"\n",
"\n",
"plt.plot(np.sort(params[idx_sort_aic[:k_best], 1]), 'o')"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8192"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(params)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8192"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(res_all0.aic)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.0797996 , 0.08011197, 0.08099224, 0.08190939, 0.08272671,\n",
" 0.08299485, 0.08720287, 0.11209125, 0.11255463, 0.11280649,\n",
" 0.11308717, 0.1133147 , 0.11356627, 0.11392479, 0.11442241,\n",
" 0.11527621, 0.11539355, 0.11541287, 0.11547739, 0.11550461,\n",
" 0.11561659, 0.11573795, 0.1157392 , 0.11595857, 0.11623736,\n",
" 0.11625448, 0.11646743, 0.11660753, 0.11668588, 0.11671932,\n",
" 0.11677643, 0.11691104, 0.11700046, 0.11710002, 0.11714394,\n",
" 0.11716887, 0.11721375, 0.11730327, 0.11744211, 0.11748089,\n",
" 0.11748219, 0.11756837, 0.1176474 , 0.11769054, 0.11780611,\n",
" 0.11786433, 0.11796896, 0.11817204, 0.11818977, 0.11832809,\n",
" 0.11838448, 0.11838829, 0.11845136, 0.11846015, 0.11851583,\n",
" 0.11854143, 0.11854223, 0.11863786, 0.11875612, 0.11897929,\n",
" 0.11899777, 0.11900155, 0.11906314, 0.11907974, 0.11910566,\n",
" 0.11933322, 0.11939592, 0.11959657, 0.11961615, 0.11963708,\n",
" 0.11966115, 0.11986583, 0.11987871, 0.12015968, 0.12041327,\n",
" 0.12054467, 0.12070152, 0.12071758, 0.12076771, 0.12103576,\n",
" 0.1210434 , 0.12143035, 0.12144246, 0.12147717, 0.12155534,\n",
" 0.1218056 , 0.12194097, 0.12200472, 0.12209712, 0.12247389,\n",
" 0.12279504, 0.12286699, 0.12288546, 0.12349604, 0.1235888 ,\n",
" 0.12443606, 0.12513315, 0.12551357, 0.12621119, 0.12681331])"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.sort(params[idx_sort_aic[:k_best], 1])"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0xf8293750c8>]"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAaYUlEQVR4nO3df6xk5X3f8ffHl4VcXMXXPzZV9gLZtYI2RUFl7WvslgqluM5CHHkRxQVSJViywj9BjVFKtahR0tI/WIu2diqhyMg/YjuJWduhdGWirCw2UivSkL1rNuAFb9hgbO5dEjaFS6L4KtyFb/+YGTw7zJl5zsw5M2ee+byk1d45d+bc55kz8z3P8z3P8xxFBGZmlq+3TLsAZmZWLwd6M7PMOdCbmWXOgd7MLHMO9GZmmTtv2gXo9a53vSt27tw57WKYmc2UY8eO/U1EbO/3u8YF+p07d7K6ujrtYpiZzRRJ3yv6nVM3ZmaZc6A3M8ucA72ZWeYc6M3MMudAb2aWucaNujGz2fDQ4+vce/gkpzc22bG0yJ17d3P9nuVpF8v6cKA3s9Ieenydux58ks2t1wBY39jkrgefBHCwbyCnbsystHsPn3wjyHdsbr3GvYdPTqlENohb9GZW2umNzaTtTu80g1v0ZlbajqXFods76Z31jU2CH6Z3Hnp8fUKltA4HejMr7c69u1nctnDOtsVtC9y5d/cbj53eaQ6nbsx6ON0wXOf9GPQ+paZ3rH4O9BVygJh9Hk2S7vo9ywPfkx1Li6z3CepFaR+rj1M3FXE+Mg9ON1QnJb1jk+FAXxEHiDw43VCd6/csc88Nl7O8tIiA5aVF7rnhcveMpsCpm4o4QOShSemGHFKBw9I7Nhlu0VckZbiZNV9T0g1OBVqVHOgr0pQAYeNpSrrBqUCrklM3FUkZbmazoQnpBqcCrUoO9BVqQoCw2dbJy0fB70dNBeaQ77fROdCbNUTvGP5eo6YCPTfAnKM3a4h+efmOca4VON9vbtHbXKkrhVHFfovy7wIe3X/NyGVzvr/56k6tOdDb3KgrhVHVfusawz9ov87dT98kUmtO3djcqCuFUdV+6xqiW7Tff/lT2z1WvwEmkVpzi34GuNVVjbpSGFXtt64hukX7HRRgcvh8zcr3ZhKpNQf6hpv2iIlZ+bKkqCuFUWXKpa4huv32e8fB432fm0PuftrfmzImseyGUzcNN80RE7lNw68rhTGrs6JzXraj6HvziYPHuerAkUZ9hifx+XGLnma3Wqc5YiK3rn1dKYxZnRV9597dbxq3PwsnqBSDvh/drXsYftzqjg+T+PzMfaBvehdvmqsp5jgsr64UxizOip7VE1SKou9Nx+bWa/ynQyf4h7OvD/zuTyo+1P35mfvUTdMnk/Tr1onWB67uLmjOXXtofYmvOnCk8uUGZsn1e5Z5dP81fPfAh3l0/zVZBHno/73ptbG5NfS73/T4kGruW/RNb7V2t7rWNzYRvBGY6u59zFrXvkwXu67lBnLUlNRmmXL0fm/K6P7uNz0+pMqyRd9pqe3a//DQVu8stFo7ra7lpcU3tT7rbF00ZcneFGUvHNe13EBumnJBfpRydL43n77pir694iLd3/1ZiA8psmvRl82pTaPVOmoLaRqti1nJPZe9oFrXcgO5acoF+XHKMahX3Kv3uz9rvdoiSS16SddKOinplKT9fX5/taRvSTor6cau7VdI+r+STkh6QtJNVRa+n7I5tUm3Wge1TIb1RHJpXdSh7EnQ72WaKhsXZXraVZdjUK+4o993f5Z6tYMMbdFLWgDuAz4ErAFHJR2KiKe6nvZ94GPAv+95+Q+AX4qIZyTtAI5JOhwRG5WUvo9RPhCTbLUWnYhSRgDk0rqoQ9nRSX4v01Q16mvc0StVlSO1Jzeo1905YU37mkUZKS36K4FTEfFsRLwKPADs635CRDwXEU8Ar/ds/4uIeKb982ngRWB7JSUv0PSWWtEHLWUEQC6tizqUnXTi9zJNVZN5xp3AlFqOcXrFndfu3P8wdxw8XtjrbsI1i7JScvTLwPNdj9eA95f9Q5KuBM4H/rLP724DbgO45JJLyu76HE1vqQ0b39ur98QwKznzSRtlTLjfy2LdLdq3LW7jR7a9hY0fbFW+BDOkte5Tjm9Kr6EoPnRmR3e2Dxr00IRrFmWlBPp+F6iL0lz9dyD9OPBl4NaIeL339xFxP3A/wMrKSql995rEJJBRLqZ2XtPvYtCgi0N19kSK6tGU4XRl1R24Z/V9Kas3YG5sbrG4bYFP3XTFyPVNmcA0LFgOO74pF2zLzI7uNehk1fThlimBfg24uOvxRcDp1D8g6UeBh4Ffj4g/LVe80dT5hR8l19j7muCHwb3MCIAqFdVj9Xsv8QfH1hs7U3hamj6Dutc4J6U6Rtr0a0n3mtQqomVmR3frNLqmNVN9HCmB/ihwqaRdwDpwM/ALKTuXdD7wP4EvRcTXRi5lg5T5EnS34nsFsCDxWvQP88s1txiL6vGVx55/U5lmoWtat1GDX929gN4UiwQv/2BrrIl1owxoGFbPlAlMnWA56ns2zgXbYT2O7kZXk1PDRYZejI2Is8DtwGHgaeCrEXFC0t2SPgIg6X2S1oCPAp+RdKL98n8DXA18TNLx9r8raqnJhKR+Cbov2hQpCvKdEQB1BtaiehSVqeld07oV1X/QUhR1X7jr3f/G5hYv/2ALGJxjHrS/UZaESK3noAlMnWA5zns2zoXjoqVG4NwL9YMu4o8zfLRuSROmIuIPgT/s2fYbXT8fpZXS6X3d7wK/O2YZGyW11ZCS8ytq0U+iG1hUj2mWqckGtfiKWsyDRprce/jkyK37QT3FQYa1yEddEqJsPQddR7vqwJFKJkaV7Q2UeW2/1E/TU3vZzYytW79cY/ciY50Px7AW8OK2Bf71e5fPyYd3tk+iG1g0+mCaZWqyYTnmfsFo3JEm/QwLyIMMOlkPWxJiUMAcpZ5F19HK9Jj7BeVxrs+N89qmzCAukk2gn9SIiNRFxga1ALu/OCs/8Y6pjHwZ1IIpKlNZOY1SSckx9wajlJEmqa37UVvxHcNO1uMsCVHFiJph++o+STWx9dz0xc8UBTnZaVlZWYnV1dVSr+nXylnctjDSBJgywemqA0f6fig7gXzUMlVZn2nJoQ5FBh333tmVqa3vQe/NqK34TiMk5cJ+ap36SSmfgO8e+PDQMqd8borK2invNBoU47x/VZF0LCJW+v0uixZ9Vd2mfi2FOw4e5xMHj7PUHtHQPWlk0Fl8nHxh07uBKcatQ5N7A6mT8lJ6AR2D3puU6z1A389o0Ymj970dZ6JhmRE1ZfZVdOzrSImNq+kTNbMI9FV1m/p9oTr9nY3NrTe2dT5MSxdue2OUQ7fOh3rUnN80uoFVB9Zx6lB113ycug16bZkLdymt3lHfszI9paL39p4bLueeGy4f+X0aVM+yAW/Y96bKVFFVJjFRcxxZBPq6FzzqZ3PrNS447y0sbluo/Cw+6dsH1pHzTK1Dv0BaZY9mnLoNe22ZsozT6k293pNi0HubMqS3zHj51IBX9kQ8iclXo6h7ZvY4srjxSFULL5UNpK9sbtWyMNYk7grfrY7bpaXUoWjMdOrFzhTj1G3U1xaNp04ZR95P0Xv56ZuuKD3fooqeVup4+ZTbE456Q5HO967IvA8H7pVFi76qblNKS6HbjqXFWs7ik+4G1pEqSqlDUSAdNpa/TAtwnLqNOkN0WA+i7PGt8vMwTm+xjmtHw/Y5bBhlFamieZBFoIe0blPZadrTWoemU5ZJdQPrShUNq8Og2blFKbGyqZg6psWXHY/eLxiOkvqp4vMwzkXDOhoEg/ZZx0lzXmUT6IsUrRqZMpGj3zoi4yzVOg0prd/UL3/VF2yLAmkn71zFzMlxAtsor236eOpxAmMdDYJB+6zrpDmPshhHXyRlpMMkx7lOWpmx7MOCeL99dU6cqcP6ev/WoCWbiy4y7tr/cN9e1qBx2nWNuumnCeOp61LH3IhB+7zj4PHSx3qeZT+Ovsi4a0zPujI51WGtorJDTzv77JW6ZHPvfjoBd5R1+8dp8ZV9bdPHU4+jjjTJoH0WjVDyhdbysg70KUE85w9NlWmEskNPy0z+KVqyuXuEy6gLbk1a7jnjugYfVJ12s3NlHejLrDGdoypzquPeAnHY9kHLI4+z4NY0OGdcjdxPmpOUdaAvWmkydf2PWVL1tPZeoww9LdpednnkcRbcstnmk2Y1spgwVaR7YkVnQtOnbrqC5xImcsySokknQGUTunonqfS7kXDHKJN/bnn/xYWTiAadNMxsuKxH3cyLaYz0GGfoadmbko8y2qPJi6KZ1WHQqBsH+gyMMuRw1pQJ3DkvkWxWZG6HV86LSS+CNkhdLekyudoclnk2q1LWOfpBmnwj37ImvQhakbpvhp2q6bNTzSZtLlv0TbwV2TiaMgxt2GqPkypfk3o4Zk0wlzn6nKepT1PRtQKg7yJlVeTM+6WK4M0TrJyjt9wNytHPZepmXrv2daerilrMC1Ll693DZIaVmuVgLlM389i1n0S6qmiCVtEkq3FPrOPeLclsXsxli74pFy/r0q/lXsddpHr1m6A26E5A455Y57VnZlbWXLbom3Lxsg5FLfe6WtW9ioZB1rE41Tz2zMxGMZeBHvJdQ2PU2/PVqa4Tq1c3NEszt4E+V6Pcnm8Scri3rtmscqDPzCi355tlufbMzKrkQJ+ZQekMB0Wz+eRAnxmnM8yslwN9htxyN7NuczmO3sxsniQFeknXSjop6ZSk/X1+f7Wkb0k6K+nGnt/9kaQNSd+oqtBmZpZuaKCXtADcB1wHXAbcIumynqd9H/gY8Pt9dnEv8IvjFdPMzEaV0qK/EjgVEc9GxKvAA8C+7idExHMR8QTweu+LI+IR4O+qKKyZmZWXEuiXgee7Hq+1t1VG0m2SViWtnjlzpspdm5nNvZRArz7bKl3EPiLuj4iViFjZvn17lbs2M5t7KYF+Dbi46/FFwOl6imNmZlVLCfRHgUsl7ZJ0PnAzcKjeYpmZWVWGBvqIOAvcDhwGnga+GhEnJN0t6SMAkt4naQ34KPAZSSc6r5f0f4CvAR+UtCZpbx0VMTOz/ubynrFmZrnxPWPNzOaYA72ZWeYc6M3MMudAb2aWOQd6M7PMOdCbmWXOgd7MLHMO9GZmmXOgNzPLnAO9mVnmHOjNzDLnQG9mljkHejOzzDnQm5llzoHezCxzDvRmZplzoDczy5wDvZlZ5hzozcwy50BvZpY5B3ozs8w50JuZZc6B3swscw70ZmaZc6A3M8ucA72ZWeYc6M3MMudAb2aWOQd6M7PMOdCbmWXOgd7MLHMO9GZmmXOgNzPLXFKgl3StpJOSTkna3+f3V0v6lqSzkm7s+d2tkp5p/7u1qoKbmVmaoYFe0gJwH3AdcBlwi6TLep72feBjwO/3vPYdwG8C7weuBH5T0tvHL7aZmaVKadFfCZyKiGcj4lXgAWBf9xMi4rmIeAJ4vee1e4FvRsRLEfEy8E3g2grKbWZmiVIC/TLwfNfjtfa2FEmvlXSbpFVJq2fOnEnctZmZpUgJ9OqzLRL3n/TaiLg/IlYiYmX79u2JuzYzsxQpgX4NuLjr8UXA6cT9j/NaMzOrQEqgPwpcKmmXpPOBm4FDifs/DPyspLe3L8L+bHubmZlNyNBAHxFngdtpBeinga9GxAlJd0v6CICk90laAz4KfEbSifZrXwL+C62TxVHg7vY2MzObEEWkptsnY2VlJVZXV6ddDDOzmSLpWESs9PudZ8aamWXOgd7MLHMO9GZmmXOgNzPLnAO9mVnmHOjNzDLnQG9mljkHejOzzDnQm5llzoHezCxzDvRmZplzoDczy5wDvZlZ5hzozcwy50BvZpY5B3ozs8w50JuZZc6B3swscw70ZmaZc6A3M8ucA72ZWeYc6M3MMudAb2aWOQd6M7PMOdCbmWXOgd7MLHMO9GZmmXOgNzPLnAO9mVnmHOjNzDLnQG9mljkHejOzzDnQm5llLinQS7pW0klJpyTt7/P7CyQdbP/+MUk729vPl/QFSU9K+nNJP1Np6c3MbKihgV7SAnAfcB1wGXCLpMt6nvZx4OWI+EngU8An29t/GSAiLgc+BPw3Se5FmJlNUErQvRI4FRHPRsSrwAPAvp7n7AO+2P7568AHJYnWieERgIh4EdgAVqoouJmZpUkJ9MvA812P19rb+j4nIs4CrwDvBP4c2CfpPEm7gPcCF/f+AUm3SVqVtHrmzJnytTAzs0IpgV59tkXicz5P68SwCnwa+BPg7JueGHF/RKxExMr27dsTimRmZqnOS3jOGue2wi8CThc8Z03SecDbgJciIoA7Ok+S9CfAM2OV2MzMSklp0R8FLpW0S9L5wM3AoZ7nHAJubf98I3AkIkLShZLeCiDpQ8DZiHiqorKbmVmCoS36iDgr6XbgMLAAfD4iTki6G1iNiEPA54AvSzoFvETrZADwY8BhSa8D68Av1lEJMzMrplZ2pTlWVlZidXV12sUwM5spko5FRN9RjR7TbmaWOQd6M7PMOdCbmWUuZXilmZmV8NDj69x7+CSnNzbZsbTInXt3c/2e3nmmk+NAb2ZWoYceX+euB59kc+s1ANY3NrnrwScBphbsnboxM6vQvYdPvhHkOza3XuPewyenVCIHejOzSp3e2Cy1fRIc6M3MKrRjabHU9klwoDczq9Cde3ezuG3hnG2L2xa4c+/uKZXIF2PNzCrVueDqUTdmZhm7fs/yVAN7L6duzMwy50BvZpY5B3ozs8w50JuZZc6B3swscw70ZmaZc6A3M8ucA72ZWeYc6M3MMueZsWZT1rSbVFh+HOjNpqiJN6mw/DjQm03RoJtUONA336z0xhzozaaoiTepsDSz1BvzxVizKWriTSosTRNvGVjEgX6GPfT4OlcdOMKu/Q9z1YEjPPT4+rSLZCU18SYVlmaWemNO3cyoWeo2WrEm3qTC0uxYWmS9T1BvYm/MgX5G+SJePpp2kwpLc+fe3ec0tqC5vTEH+hk1S91GsxzNUm/MgX5GzVK30SxX4/TGJjk00xdjZ5Qv4pnNrs41tvWNTYIfXmOra0CFA/2Mun7PMvfccDnLS4sIWF5a5J4bLm9kt9HMzjXpoZlO3cwwX8Qzm02TvsaWFOglXQv8FrAAfDYiDvT8/gLgS8B7gf8H3BQRz0naBnwWeE/7b30pIu6psPxmZo3WLxc/6WtsQwO9pAXgPuBDwBpwVNKhiHiq62kfB16OiJ+UdDPwSeAm4KPABRFxuaQLgackfSUinqu6ImZmTdNvvssdB48TgIDoem6d19hScvRXAqci4tmIeBV4ANjX85x9wBfbP38d+KCkTj3eKuk8YBF4FfjbSkpuZtZw/XLx0fW/2j/XfY0tJXWzDDzf9XgNeH/RcyLirKRXgHfSCvr7gBeAC4E7IuKl3j8g6TbgNoBLLrmkZBXMzJppWM49aAX5R/dfU2s5Ulr06rMtEp9zJfAasAPYBfyapHe/6YkR90fESkSsbN++PaFIZmbNl5Jzn8Qkx5RAvwZc3PX4IuB00XPaaZq3AS8BvwD8UURsRcSLwKPAyriFNjObBf3mu/SaxCTHlEB/FLhU0i5J5wM3A4d6nnMIuLX9843AkYgI4PvANWp5K/AB4DvVFN3MrNm657vAm1Mfk5rkODRH38653w4cpjW88vMRcULS3cBqRBwCPgd8WdIpWi35m9svvw/4AvBtWnX8QkQ8UUM9zMwaqXu+y7TuSKVWw7s5VlZWYnV1ddrFMDObKZKORUTf1LiXQDAzy5wDvZlZ5hzozcwy50BvZpY5B3ozs8w1btSNpDPA98bYxbuAv6moOLNiHusM81nveawzzGe9y9b5JyKi79ICjQv045K0WjTEKFfzWGeYz3rPY51hPutdZZ2dujEzy5wDvZlZ5nIM9PdPuwBTMI91hvms9zzWGeaz3pXVObscvZmZnSvHFr2ZmXVxoDczy1w2gV7StZJOSjolaf+0y1MXSRdL+mNJT0s6IelX29vfIembkp5p///2aZe1apIWJD0u6Rvtx7skPdau88H2/RKyImlJ0tclfad9zP9Z7sda0h3tz/a3JX1F0o/keKwlfV7Si5K+3bWt77Ft39Pjf7Tj2xOS3lPmb2UR6CUt0Fr7/jrgMuAWSZdNt1S1OQv8WkT8E1o3cvmVdl33A49ExKXAI+3HuflV4Omux58EPtWu88vAx6dSqnr9Fq27tP0U8E9p1T/bYy1pGfh3wEpE/DSte2DcTJ7H+neAa3u2FR3b64BL2/9uA367zB/KItDTujftqYh4NiJeBR6gdVPy7ETECxHxrfbPf0fri79Mq75fbD/ti8D10ylhPSRdBHwY+Gz7sYBraN2AHvKs848CV9O6sQ8R8WpEbJD5saZ1Q6TF9m1JLwReIMNjHRH/m9aNmroVHdt9wJei5U+BJUk/nvq3cgn0y8DzXY/X2tuyJmknsAd4DPjHEfECtE4GwI9Nr2S1+DTwH4DX24/fCWxExNn24xyP+buBM8AX2imrz7ZvyZntsY6IdeC/0roN6QvAK8Ax8j/WHUXHdqwYl0ug770VI0DW40Yl/SPgD4BPRMTfTrs8dZL088CLEXGse3Ofp+Z2zM8D3gP8dkTsAf6ejNI0/bRz0vuAXcAO4K200ha9cjvWw4z1ec8l0K8BF3c9vgg4PaWy1E7SNlpB/vci4sH25r/udOXa/784rfLV4CrgI5Keo5WWu4ZWC3+p3b2HPI/5GrAWEY+1H3+dVuDP+Vj/K+C7EXEmIraAB4F/Tv7HuqPo2I4V43IJ9EeBS9tX5s+ndfHm0JTLVIt2bvpzwNMR8d+7fnUIuLX9863A/5p02eoSEXdFxEURsZPWsT0SEf8W+GPgxvbTsqozQET8FfC8pN3tTR8EniLjY00rZfMBSRe2P+udOmd9rLsUHdtDwC+1R998AHilk+JJEhFZ/AN+DvgL4C+B/zjt8tRYz39Bq8v2BHC8/e/naOWsHwGeaf//jmmXtab6/wzwjfbP7wb+DDgFfA24YNrlq6G+VwCr7eP9EPD23I818J+B7wDfBr4MXJDjsQa+Qus6xBatFvvHi44trdTNfe349iStUUnJf8tLIJiZZS6X1I2ZmRVwoDczy5wDvZlZ5hzozcwy50BvZpY5B3ozs8w50JuZZe7/A3CoYgeGCtWBAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(params[idx_sort_aic[:k_best], 1], 'o')"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0xf8293cecc8>]"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD4CAYAAAAHHSreAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAem0lEQVR4nO3df7Bc5X3f8feXq4u5hjZXCiIDF2TJsUIKJpHCIuzSMhTXSDg2UgiuwTQRmaTqH2USOzEdKdNGRrhGCXFxOmVIsC2XdGrAg7Gi2Ew0sjHTxpMyukIEkLEGDZHRvaIgR1zHsa6Nfnz7x56Fo91zds/Zn8855/OauSPt2T1nz9mz+3yf34+5OyIiInFnjPoEREQkPAoOIiLSQsFBRERaKDiIiEgLBQcREWmxYNQn0A/nnnuuL126dNSnISJSKHv27Pm+uy9Oeq4UwWHp0qVMT0+P+jRERArFzL6X9pyqlUREpIWCg4iItFBwEBGRFgoOIiLSQsFBRERalKK3kohIqLbvneWenfs5PDfPBZMT3LH6YtatnBr1aXWk4CAiMiDb986y6bHnmD9+EoDZuXk2PfYcQPABQtVKIiIDcs/O/W8Ghob54ye5Z+f+EZ1RdgoOIiIDcnhuPtf2kCg4iIgMyAWTE7m2h0TBQURkQO5YfTET42OnbZsYH+OO1ReP6IyyU4O0iMiANBqd1VtJREROs27lVCGCQTNVK4mISAsFBxERaaHgICIiLRQcRESkhYKDiIi0UHAQEZEWCg4iItJCwUFERFooOIiISAsFBxERaaHgICIiLRQcRESkhYKDiIi0UHAQEZEWmYKDma0xs/1mdsDMNiY8f7WZPW1mJ8zspqbn1pvZi9Hf+tj2J6NjPhP9nRdtv83MjsS2/1avFykiIvl0XM/BzMaA+4D3AzPAbjPb4e7fib3sZeA24BNN+y4CNgM1wIE90b6vRy+51d2nE972EXe/Pe/FiIhIf2QpOawCDrj7S+7+BvAwsDb+Anc/6O7PAqea9l0N7HL3o1FA2AWs6cN5i4jIAGUJDlPAodjjmWhbFp32/WJUdfSfzcxi23/VzJ41s0fN7KKkA5vZBjObNrPpI0eOZDwdERHJIktwsIRtnvH47fa91d0vA/5l9Pdr0fa/BJa6+y8A3wAeTDqwuz/g7jV3ry1evDjj6YiISBZZgsMMEM+9Xwgcznj81H3dfTb694fAl6hXX+Huf+/uP4le/zng8ozvJSIifZIlOOwGlpvZMjM7E7gZ2JHx+DuB68xsoZktBK4DdprZAjM7F8DMxoEPAs9Hj8+P7X8D8ELG9xIRkT7p2FvJ3U+Y2e3UE/oxYJu77zOzLcC0u+8wsyuArwILgQ+Z2Z3ufqm7HzWzu6gHGIAt0bazqQeJ8eiY36BeSgD4bTO7ATgBHKXeC0pERIbI3LM2H4SrVqv59HRSj1gRGYbte2e5Z+d+Ds/Nc8HkBHesvph1K7P2W5FRMbM97l5Leq5jyUFEpJ3te2fZ9NhzzB8/CcDs3DybHnsOQAGiwDR9hoj05J6d+98MDA3zx09yz879Izoj6QcFBxHpyeG5+VzbpRgUHESkJxdMTuTaLsWg4CAiPblj9cVMjI+dtm1ifIw7Vl88ojOSflCDtIj0pNHorN5K5aLgICI9W7dySsGgZFStJCIiLRQcRESkhYKDiIi0UHAQEZEWCg4iItJCvZVERApo0JMdKjiIiBTMMCY7VLWSiEjBDGOyQwUHEZGCGcZkhwoOIiIFM4zJDhUcREQKZhiTHapBWqRLWhpTRmUYkx0qOIh0QUtjyqgNerJDVSuJdEFLY0rZKTiIdEFLY0rZKTiIdEFLY0rZKTiIdKFKS2Nu3zvLVVufYNnGr3PV1ifYvnd21KckQ6AGaZEuVGVpTDW8V5eCg0iXqrA0ZruG97Jfe9WpWklEUqnhvboUHEQklRreq0vBQURSVanhfdhCb+hXm4OIpKpKw/uwFaGhX8FBRNqqQsP7sBWhoT9TtZKZrTGz/WZ2wMw2Jjx/tZk9bWYnzOympufWm9mL0d/62PYno2M+E/2dF21/m5k9Er3XU2a2tLdLFOks9CK+lEsRGvo7BgczGwPuA64HLgFuMbNLml72MnAb8KWmfRcBm4ErgVXAZjNbGHvJre6+Ivp7Ldr2m8Dr7v4u4F7gD3NflUgOjSL+7Nw8zltFfAUIGZQiNPRnKTmsAg64+0vu/gbwMLA2/gJ3P+juzwKnmvZdDexy96Pu/jqwC1jT4f3WAg9G/38UeJ+ZWYbzFOmKJtGTYStCQ3+WNocp4FDs8Qz1kkAWSfvGK9S+aGYnga8An3J3j+/j7ifM7AfATwPfjx/YzDYAGwCWLFmS8XREWoVaxNd6EeVVhIb+LMEhKdfuGY/fbt9b3X3WzP4J9eDwa8CfZ30/d38AeACgVqtlPR+RFhdMTjCbEAhGWcQvQm8W6U3oDf1ZqpVmgItijy8EDmc8fuq+7j4b/ftD6m0Vq5r3MbMFwE8BRzO+n0hunYr4o2isrnJVlzoHhCFLyWE3sNzMlgGzwM3ARzMefyfw6Vgj9HXApijRn3T375vZOPBB4BvRa3YA64G/AW4Cnoiqm0QGol0Rv585+DzVRKFWdQ2aSkzh6Bgconr/26kn9GPANnffZ2ZbgGl332FmVwBfBRYCHzKzO939Unc/amZ3UQ8wAFuibWcDO6PAMEY9MHwues0XgP9pZgeolxhu7uP1iiRKK+L3qz963kRvkFVdIbdlFKH/f1VkGgTn7o8Djzdt+4PY/3dTrzJK2ncbsK1p24+Ay1Ne/2Pgw1nOS2TQ+pWDz5vo3bH64tOCCfSnN0voOfOqlphCpLmVRNroV3/0vIneupVT3H3jZUxNTmDA1OQEd994Wc8JeOhtGUXo/18Vmj5DpI1+5eC7qSYaRG+W0HPmgyoxSX4qOYi00a8cfCiDnkLPmQ+qxCT5WRk6AtVqNZ+enh71aYi0FUJDcHObA9SDlBLgajKzPe5eS3pO1UoiQxLCoKcijMyVMCg4iFRMCEGqjEIoGfaTgoOUXtl+tBKe0LsId0PBQUqtaD9aBbJiadyvpJ5oRR+8p95KUjh55t4JvV9/nNaVKJb4/UozOzdf2PmhFBykUPImoKH3648rUiCT5PuVpKhBXtVKUihZp6FoFPfTOmqH0q8/rkiBrEzSqvI6VfHluS9FrGJScJBCyZKAJvXljwt1xG2I60qUXVqb1PT3jvKVPbNt26rS7leaogV5VStJoWQZ4duuuD/sEbd52kdCGUVdJWkl0YeeOtSxii/tfi18+3jie7UL8iGuYaGSgxRKlrl30nJoBnx747U9vX+e3kR5e0ppgNrwpX1XTqbMHBF/fdr9AnLNDxVqjzoFBymULAnooKpn8v6Iu1mbQAPUhivtuzJmlhggmr9D7e5X1iAf6hoWCg7SlVH2x++UgA5qZs+8P2I1MLcKbRxH2nflVy+fOq3NobE963coT5AP9Xui4CC5hVoMbhhU9UzeH3EZG5h7SdxD/N60+67U3rFoKIEs1O+JgoPkFmoxOG4Q1TN5f8RlW5ug18Q91O9N2ndlWFV8oX5P1FtJcgu1GDxoeXsTlW1tgl4H6XXzvQmxF0+/hfo9UclBcgu1GDzo+uxuqqvK1MDca6Yg7/cmpGqoYXy3QvueKDhIbiEWg4eVkIT4Ix6WXjMFeb83oVRDhRSkhknVSpJbiMVgzUs0eL0O0sv7vQml+rKq3y2VHKQroeWgQ0lIyqwfvcDyfG9Cqb5M+w41ZlwddXfcQVFwkFIIJSEpu2FmCkKpvmw3h1KZq5gUHCQo3Tb8ZU1I4sf/qYlxzGDu2PEgBmTJ6UKZTiTpuxUXQnfcQVBwkGD00vCXJSFpPv7c/PE3nytzDrDIeimpdJPRaLdP2opvUM7qSwUHCUavvVM6JSSdFmcpaw6wirrJaHTaZ93KKa7a+kRlqi/VW0mCMehG5SzHKWMOsIq66WGUZZ8qTauu4CDByLJWwyCOP4j3ktHqJqORZZ8Qu3EPiqqVJBiD7p3SqWGxrDnAKulledisPd5C68Y9KAoO0rN+TS0w6N4pzccPpbdSaNNYF1XW5WHTPu9Qus42G9X3wzxlxaPTXmS2BvgTYAz4vLtvbXr+auCzwC8AN7v7o7Hn1gP/KXr4KXd/sGnfHcA73f3d0eNPAv8OOBK95Pfd/fF251er1Xx6errjdYSiTIlB0g9yYnysL0XtMn1OaQb5+VVNWmMxwGSUEXj92HEMTitZxD/v0L5zg/5+mNked68lPdex5GBmY8B9wPuBGWC3me1w9+/EXvYycBvwiaZ9FwGbgRr1+7En2vf16PkbgX9MeNt73f2PO51bERV5npakH04/5r9JOi5Q2M8pj1DmDyqDdu0JPzlx6s3PuTk7HP+8Q6syGuX3I0u10irggLu/BGBmDwNrgTeDg7sfjJ471bTvamCXux+Nnt8FrAEeMrNzgN8FNgBf7u0yiqOoiUFaUEsrwmft9ZN23LPGzyjk55RXGaf9GEbuO+k92i352a4LM4T1ecevLa1eZxhTd2QJDlPAodjjGeDKjMdP2rdxJXcBnwGOJex3u5n9OjAN/F6jpBFnZhuoBxaWLFmS8XRGr6iJQVpQy7rWbt7j9hp0RilP4hjqtB/dJvDDKBmnvUfa0p6dAgOM/vNu6NRuEjfo0nSWrqyWsK1zQ0Wbfc1sBfAud/9qwvP3Az8LrABeoR5AWg/i/oC719y9tnjx4oynM3qD7q45KGmJ8kn3ln7fxls5m06Ls+RN7EP/nBo/7tko19f4Aad9DiH2m897DXHDmME07T2+9d0jid1Mpzp8Z0b9ecd1GqjZbJCzw2YJDjPARbHHFwKHMx4/bd/3Apeb2UHgr4GfM7MnAdz9VXc/6e6ngM9Rr9YqjRATgyzSEuXmH2C8sS9LopJ23MmJ8UJ+TnkTxxD7zfeSwA+jZNzuPdatnOLbG6/l77b+Mt/eeO2bvZCSMjAQxucd183nNKjSdJZqpd3AcjNbBswCNwMfzXj8ncCnzWxh9Pg6YFPUBnE/gJktBb7m7tdEj89391ei1/8K8HzG9yqEUCYTy6tdN792Uwt0aidIO+4nb7gUKN7n1E3iGFojaC8J/DCqyfK+R5F+c2nX1sh8DbMKsmNwcPcTZnY79YR+DNjm7vvMbAsw7e47zOwK4KvAQuBDZnanu1/q7kfN7C7qAQZgS6Nxuo0/iqqdHDgI/PvuLi1coSUGWWT5gXWbMLY7btE+p1DbEPLo5RqSgv34GcaxN06wbOPX+5IwdzMeIZTfXKe2nE7XNsxxGJnGOYSuaOMcyiqtn/nU5ATf3njtCM5o+MowbqHXa2ieFv1Hb5zg+Mm30pl+fB6hjUfIIuvn2u7a+n3d7cY5KDhIW3m+jGVIGPuhiAlXs35dgzIMbwnxs+hpEJxUV95uiaOo2w0xIQ6lCqMX/bqGonbdHoSifRYKDgL0b/TzoBLGYY+iDjHopAn5XMvQBtMvRfssNGW3pPZrD2XVq7Tzu/Mv9w2kT30v/fyHLfRzLWrX7UEo2meh4FAB2/fOctXWJ1i28euJA9PajX5OMuycTtr5vX7seOLrew1ewxjI1S9p5/qxR57JNAhx0EIcxzEqRfssVK00JKMq+mdpN+g0+nnUUxgPexR1keqG251TKJMVlqENpl+K9Fmo5DAEoyz6Z8kFZxn9PMqczrBHUQ9jipNOpbmsOp1TqCUeCZ+CwxCMspoiSy64XV1o0nQEw5Z2fp+84dKBBK9B1w33M7OQdK7NQizxSPhUrTQEo6ymyNJDIvTpBYY9inrQn0c/p22Pn2taB4JQe8NI2BQchmCUXdiyTjUw7LrQvG0wwz6/Qb5fvzMLjXNNG4QYam8YCZuCwxCMcm3akEoFjYAwOzefOHtr/HzLbFCZhZDutRSfps8YkpAHKg1Kpzl2mlVlSgVNMyKh0PQZAShSF7Z+aE4A5+aTxyTEVaXhVDn8/qpixmsYFBxkIPKuaAXVajitWmahk5CXJa0qBQcZiLylgKI3nCr32r1eEvh+9vyS02mcgwxEllJAqEs15hX6/EahC31Z0qpSyUFa9CMXnLYi2DlnLWDu2PFS5a6Ve+1OvPdaklCWJa0qBQc5Tb/qcKvU6Krca35JPbaaNSfwSZmWUXYTLzsFBzlNv0fvljEYNAs19xpyO0inDgvNCXxapuXuGy/j7hsvC/Y6i0zBoaSyJAxJr1EuOL9R5F473d9h9OLpJfi0+z5NJRyrXaZlVHN+lZ2CQwllSRjSXjP59vHEdRJGnQsO2bCr0LLc33brPNyzc/+b5zeqLqRppa20gZDKtAyfgkMJder9kdYIOH/8JG9bcEYQazgUzTCr0LJU/WVZ52H6e0f5yp7Z0xL4jz/yDB975JnE3Hvec2gnb2kr1Kq7MlNX1hJKSxg6Lf8J8IP540Gs4SCni6//kKV3T5Z1Hh566lBLAt8831Vad9xec/J5V0Ur2hKbZaCSQwml5bLGzDqOWr5gcqIyDclFkaVnD5weEJJy5s1OdphXrV1JoB85+Tzfsyr1fguFgkOXQu4JklZk75S4KCcWpixTkTTfuyzrPIyZdQwQaSWBUTTCK9MyXKpW6kLoI2LTiuxTbXJ1qj4KV7uqmnZVMo1V/D77kRWJVTK3XHlRx1Xk0koCeauFpHhUcuhCEUbEpuWyNFV08eTt2dOsXZVM7R2LEtfYgM4lAeXky03BoQtZG+NCq3pSvW0x9aMKJy0hj28P7fsqo6Xg0IUsjXGjmEo4y49bub3iGVZQ13dD4hQcupAlJzfsqifNa19uSrhl2BQccmhe9vKs8TNSZxgd9ojOIrSDiEhxZOqtZGZrzGy/mR0ws40Jz19tZk+b2Qkzu6npufVm9mL0tz5h3x1m9nzs8SIz2xW9fpeZLezmwvqtuYfS3Pxxfnz8FPd+ZEXi3C5pvTwGNaJT0wuISD91DA5mNgbcB1wPXALcYmaXNL3sZeA24EtN+y4CNgNXAquAzfHE3sxuBP6x6VgbgW+6+3Lgm9Hjkcu7IMmwR3QOOxiJSLllKTmsAg64+0vu/gbwMLA2/gJ3P+juzwKnmvZdDexy96Pu/jqwC1gDYGbnAL8LfKppn7XAg9H/HwTW5biegcmbMx92P3BNLyAi/ZSlzWEKOBR7PEO9JJBF0r6N1PEu4DPAsaZ9fsbdXwFw91fM7LykA5vZBmADwJIlSzKeTve6mS5gmI2I6qYqIv2UJThYwrb2Y+477GtmK4B3ufvHzWxpxmOdfhD3B4AHAGq1Wtbz6VoRVpxSjxYR6ZcswWEGuCj2+ELgcMbjzwDXNO37JPBe4HIzOxidw3lm9qS7XwO8ambnR6WG84HXMr7XQGXNmWsgkYiUQZbgsBtYbmbLgFngZuCjGY+/E/h0rBH6OmCTux8F7geISg5fiwIDwA5gPbA1+vcvMr7XwHXKmWusgYiURccGaXc/AdxOPaF/Afiyu+8zsy1mdgOAmV1hZjPAh4E/M7N90b5Hqbct7I7+tkTb2tkKvN/MXgTeHz0uhLw9mkREQpVpEJy7Pw483rTtD2L/3029yihp323AtjbHPgi8O/b474H3ZTmv0GisgYiUhUZI90GjnSGtVVxjDUSkaBQcetRpla7QejSJiGSh4NCjdqt0dVqkXUQkVAoOPUprTzDItBCLur6KSIi0TGiPepnTKPTlRkWkuhQcetTLnEbq+ioioVK1Uo96mdNIXV9FJFQKDn3Q7ZxG3UzmJyIyDKpW6mD73lmu2voEyzZ+nau2PtHX9gBNsy0ioVLJoUnzUqA/euMEx0/Wh7f1e66kbqqk1LtJRIbB3Ac+2/XA1Wo1n56e7vk4nQa0NUxNTmTqptpvSec3MT420EWERKS8zGyPu9eSnlO1Uky7AW1xo2owTuvd9LFHnul7lZeIVJuCQ0zWRH9UDcbtzk9jJESknxQcYrIk+qNoMG40ineqANQYCRHpFzVIxyQtBTp+hnHOWQuYO3Z8qA3AjYbn2bl5jOzrsmqMhIj0g4JDTC8D2vqpueE5T5cBjZEQkX5QcGjS7YC2fsraMD4xPtbSc0ljJESkH9TmEKAsVUNTkxPcfeNlTE1OYLHHow5sIlIOKjkEKG1ajYZGCSGEUo6IlJOCQ0DaNUI3HmsBIREZBgWHQCQ1QisgiMioKDgEIqkRuhEYRjFVh4hUmxqkA6G1HUQkJAoOgehluVERkX5TcAiE1nYQkZCozSEQoYzOFhEBBYegaNyCiIRC1UoiItJCwUFERFooOIiISAsFBxERaZEpOJjZGjPbb2YHzGxjwvNXm9nTZnbCzG5qem69mb0Y/a2Pbf8rM/tbM9tnZn9qZmPR9k+a2ayZPRP9faDXixQRkXw6Boco0b4PuB64BLjFzC5petnLwG3Al5r2XQRsBq4EVgGbzWxh9PS/cfdfBN4NLAY+HNv1XndfEf09nvuqRESkJ1m6sq4CDrj7SwBm9jCwFvhO4wXufjB67lTTvquBXe5+NHp+F7AGeMjd/yF2DmeSb8GznjVmQNWYAhGRVlmqlaaAQ7HHM9G2LNrua2Y7gdeAHwKPxl53u5k9a2bbYiWN05jZBjObNrPpI0eOZDydusYMqLNz8zgwOzfPpseeY/ve2VzHEREpqyzBwRK2Zc3lt93X3VcD5wNvAxpTj94P/CywAngF+EzSgd39AXevuXtt8eLFGU+nLmkG1PnjJ7ln5/5cxxERKasswWEGuCj2+ELgcMbjd9zX3X8M7KBeVYW7v+ruJ939FPA56tVafaUZUEVE2ssSHHYDy81smZmdCdxMPTHPYidwnZktjKqHrgN2mtk5ZnY+gJktAD4AfDd6fH5s/18Bns/4XplpBlQRkfY6Bgd3PwHcTj2hfwH4srvvM7MtZnYDgJldYWYz1Hsc/ZmZ7Yv2PQrcRT3A7Aa2RNvOBnaY2bPA31Jvd/jT6C3/yMyei577V8DH+3e5dZoBVUSkPXMfaiehgajVaj49PZ1rH/VWEpGqM7M97l5Leq6ys7JqBlQRkXSaPkNERFooOIiISAsFBxERaaHgICIiLRQcRESkhYKDiIi0qGxX1jiNeRAROV3lg0NjhtbGRHyNGVoBBQgRqazKVytphlYRkVaVDw6aoVVEpFXlg4NmaBURaVX54KAZWkVEWlW+QbrR6KzeSiIib6l8cADN0Coi0qzy1UoiItJKwUFERFooOIiISAsFBxERaaHgICIiLczdR30OPTOzI8D3utz9XOD7fTydoqjidVfxmqGa113Fa4b81/0Od1+c9EQpgkMvzGza3WujPo9hq+J1V/GaoZrXXcVrhv5et6qVRESkhYKDiIi0UHCAB0Z9AiNSxeuu4jVDNa+7itcMfbzuyrc5iIhIK5UcRESkhYKDiIi0qHRwMLM1ZrbfzA6Y2cZRn88gmNlFZvYtM3vBzPaZ2e9E2xeZ2S4zezH6d+Goz7XfzGzMzPaa2deix8vM7Knomh8xszNHfY79ZmaTZvaomX03uufvrci9/nj0/X7ezB4ys7PKdr/NbJuZvWZmz8e2Jd5bq/tvUdr2rJn9Ut73q2xwMLMx4D7geuAS4BYzu2S0ZzUQJ4Dfc/d/BrwH+A/RdW4Evunuy4FvRo/L5neAF2KP/xC4N7rm14HfHMlZDdafAH/l7j8P/CL16y/1vTazKeC3gZq7vxsYA26mfPf7fwBrmral3dvrgeXR3wbg/rxvVtngAKwCDrj7S+7+BvAwsHbE59R37v6Kuz8d/f+H1BOLKerX+mD0sgeBdaM5w8EwswuBXwY+Hz024Frg0eglZbzmfwpcDXwBwN3fcPc5Sn6vIwuACTNbALwdeIWS3W93/9/A0abNafd2LfDnXvd/gUkzOz/P+1U5OEwBh2KPZ6JtpWVmS4GVwFPAz7j7K1APIMB5ozuzgfgs8B+BU9Hjnwbm3P1E9LiM9/udwBHgi1F12ufN7GxKfq/dfRb4Y+Bl6kHhB8Aeyn+/If3e9py+VTk4WMK20vbrNbNzgK8AH3P3fxj1+QySmX0QeM3d98Q3J7y0bPd7AfBLwP3uvhL4ESWrQkoS1bOvBZYBFwBnU69WaVa2+91Oz9/3KgeHGeCi2OMLgcMjOpeBMrNx6oHhf7n7Y9HmVxvFzOjf10Z1fgNwFXCDmR2kXl14LfWSxGRU7QDlvN8zwIy7PxU9fpR6sCjzvQb418DfufsRdz8OPAb8c8p/vyH93vacvlU5OOwGlkc9Gs6k3oC1Y8Tn1HdRXfsXgBfc/b/GntoBrI/+vx74i2Gf26C4+yZ3v9Ddl1K/r0+4+63At4CbopeV6poB3P3/AYfM7OJo0/uA71Diex15GXiPmb09+r43rrvU9zuSdm93AL8e9Vp6D/CDRvVTVpUeIW1mH6CeoxwDtrn7fxnxKfWdmf0L4P8Az/FW/fvvU293+DKwhPqP68Pu3tzYVXhmdg3wCXf/oJm9k3pJYhGwF/i37v6TUZ5fv5nZCuqN8GcCLwG/QT0TWOp7bWZ3Ah+h3jtvL/Bb1OvYS3O/zewh4Brq03K/CmwGtpNwb6Mg+d+p9246BvyGu0/ner8qBwcREUlW5WolERFJoeAgIiItFBxERKSFgoOIiLRQcBARkRYKDiIi0kLBQUREWvx/F8eFksl93UcAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(bse[idx_sort_aic[:k_best], 1], 'o')"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0xf829433448>]"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYlElEQVR4nO3df6zddX3H8deby0WuOL1qq7O31GJCUCZK8UbZuiiiCz800uESJWa6BdN/dNPGsZW4TOey0AUznQnDVEXELAhT1jXo1i2UhQQn43ZFQKHaiUhv0V4jF6fcjLa+98c5B849/f4839+f83wkTXt+9JzP93y/39f5fD+/jrm7AADdd1LTBQAAlINAB4BAEOgAEAgCHQACQaADQCBObuqN16xZ4xs3bmzq7QGgk/bt2/dTd18b9Vhjgb5x40YtLCw09fYA0Elm9mjcYzS5AEAgCHQACASBDgCBINABIBAEOgAEorFRLgC6Ydf+RV2754AOL69o3eyMrrroLG3ZNNd0sRCBQAcQa9f+RV192wNaOXpckrS4vKKrb3tAkgj1FqLJBUCsa/cceCbMB1aOHte1ew40VCIkoYY+Bi5Bw8B+THd4eSXX/WgWgZ4Tl6BhYD9ms252RosR4b1udmbVbb4c2yG1ycXMbjCzI2b2YMzj7zGz+/t/vmlmry2/mO3BJWgY2I/ZXHXRWZqZnlp138z0lK666Kxnbg++HBeXV+R69stx1/7FmkuLLG3oN0q6OOHxRyS9yd1fI+mvJO0soVytxSVoGNiP2WzZNKdrLj9Hc7MzMklzszO65vJzVtW++XJsj9QmF3e/y8w2Jjz+zaGb35K0vnix2ivrJSjajf2Y3ZZNc4nNJ3w5tkfZo1yulPQvcQ+a2VYzWzCzhaWlpZLfuh5ZLkHRfm3aj7v2L2rzjr06Y/vXtXnH3s41VcR9CfLlWL/SAt3M3qxeoP9Z3HPcfae7z7v7/Nq1kcv5tl6WS1C0X1v2Ywjtz236cpx05u7pT+o1udzu7q+Oefw1kv5J0iXu/r0sbzw/P++sh45Jt3nH3simn7nZGd29/cIGSjQeRrnUx8z2uft81GOFhy2a2QZJt0n6/axhDqAnlPbntHZ21CM10M3sZkkXSFpjZockfUzStCS5+2cl/YWkF0v6ezOTpGNx3x4AVqNzFmXKMsrlipTH3y/p/aWVCJgAgyaKxeUVmaThhs9x259p9gAzRYGajc5SdemZUJ8bM4iZ+QqJQAcyK6sGHDURZxDm43aEJk3uIdAnB4EOZFBmDbiKjtBQOldRDMvntkjXJ5iErMzp7VVMxGFyDyQCvTVCmGDSJmV/OZZZA65iIk7Sa1JRmBwEekuwwFF5qvhyLLMGXMUs1bjXlERFYYJ0tg19nA6qNg/raroNtM2fTV5VdBBeddFZq9rQpWK16iom4kS95uYde4PvLO3SsVt1WTsZ6ON0ULV9WFeTE0za/tnkVcWX4+Bz6EpwDDRdUahal47dOsraySaXcZon2t6k0eQCR23/bPKqqoNwy6Y53b39Qj2y4226e/uFrQuMKKF3lmY9dtvQj1DHedbJQB+n1tH2mspoG+jszLROnT5J2265r/IDsO2fTV50ED4r9JUQ447RxeWVZ/Zv1j6Vqo+NOs6zTja5jNM8UXeTxjhtZYM20LovI0NbTySueURSZy7Py9LVpqKs4o5d6dn9e+r0San9CHWcc3WcZ50M9HE6qMru1EqSdnCkhX3ds/7q/GzKkOXLcriDcHjdlFGhdRBGydoB24bOxbxliDp2h60cPR772HDNuI5zro7zrJOBPk6to86aSlpbWVpNoO4mkC7V4vLWpEafH6WrTUtlakPn4jhlGD5242rqcYZrxnWcc3WcZ5l+4KIKIf/AxRnbv66oT9UUf9k1vI5HKD96UIW8n03c87P830nShmOuaBmy7OuBmempVWP/27D9WSX9wEUnO0XbatCpEvcVuW52JlNNIPSOrCLy1qTSalh8rj1V1VDzdDQWLUPUeRMlaiJXKOdcJ5tc2ijt0n76JNNTTx9LDPuBLjWB1C1vx1JSp9m4S9WGqMwOu7i13tOaUIqWIUvzi0mRNe5QzjkCvaCkDreB2Zlp/fLpY3riqaORj0fVBPhJr2h5O5bins8Pe69WVodd1Frvw5I6Gssow+C8iWtCSfpyCOGcI9ATpPW4Z+lwM0mnPedkLa9Ehzm1xHzy1qRCqXnFKWtkSlmfU9SAgFFxTSh5ypC23Vm/HNowsqdMdIrGiArrrB0pw+b67eZxnaSP7HhbiaU+UdwBG9qBPEnSfr6uyquPtOMmbkDAsKIdjVnOzaSyNvn5lSGpU7RTNfQ6QyjLuNSsHW5xTTJVT9yJGwa28OjP9LV9ixM1wSaPNn/ZFWnSGPz/cbcty7DCpD4LafWM3XHLkXXMeFQTStHPr+06M8ql7vXCs/S4JwXycE96Uz3ocQf+zfc8FtTaLWVq07r0USNEijRpFN22LGuRRB3r1v+7rCV9i4yGKfL5dUFnAr3uBaSyLGoUF9Sffte5qxZvqmL96yziDszjMc1sXT6Qy9KWhcriwjfLOOu4Yzdu2z5y67cLDSscXjcl6lj/1LvO1Q+HFjQr+hkXWXAsyzHe1SUvpA41uZQ9TraMTpU8nThxPehVXt7HXf5OmUWGepcP5LKMc5wl7cNx929c6KVJuvJL+4Ifd1hh1P9N2sYyxpuPOxqmjiahJnWmhl7mMqBRtZ9tt9ynjUO1lKyrHxZZUrXqy/u4K4gr3nB6EJMoqhB3PLkUWYNN2ofj7N9BM0ueaeyjTRpxx2CWcyWpppw2cSdrLTvPuRzV7FTkireOJqEmdWaUS9ae7SzSTphBz/dgSKGk0t47SznKnG5cxyiXrtZmoqQNRc060inuKkiK379ZhsFGvVbWzzvr6yeNvkqbd5Fl5FaeUSpVnHdpx2vSeTkY5NDksZ40yqUzgS6VFxxZhlYNzExP6dTpkyInBRUN3qQ1X6oezliWqk66JqWF1vB+z3MsDcTt33Fq5nmPk+Fz6KSYL50sx3XRykiWcznp86hy/kbSPp2Znmr8WA9m2GJZM7nS2tGGZV1+s8xydKktu+6lfpOUOclmy6a52BN7dKRTkVX+4l63yOskGV1WeNy26KKzOrOcy0mfR5VDbZP6ntpyrMfpTBt6mbIu4pOmaPCGsCBQ0Q6usn4lpmh/RFQ5xh3plCRp/8a93+zMdCXHSVJbdNp+qWPkVtr5VdXoo7jzsgujwzpVQ88qraY2uojP6GyxUbMz0/q/Y78qfWH6uqelV9HWXeQqo8w1uItcKcSV452vm1s1AUtKH+kU14whpTcTxNV6P/6O31j1HmUeJ1km38Ttl7xXzGX/eIVUTZjGnZdNTRDMI5hAz7vC2+ilZ9JU4LpPqCpU9QMGRdbMKLO5puzJJitHj+vOh5d0zeXn5OpULtKnkPYFX2V4Zvk8ijQtVPXjFVWFadx52fZf9goi0ItO540K96InVNtU1dad5Soj7mQus2+iyJVC0pfBuDXYccO0jC/4ol/eVayNnjaZKOmcS/qirDNM676iHkcQgV7mdN66asxlS6uR5TlJ89bu0j6zuJM5aYJTGZfnRSeb5J1xOfhybPoYKvrlXUVnfdIs06zNO1LzYdr0vk0TRKCHPp03TZFFk0Y/l6jX2nbLffrwLfdpdmZaZtLyU0dznVBJMxSjmrje/Mq1hS7P857web8M6v7N17yanIkZp4yRI20P0zYIYpRLWli3rZ2rbOMumhT1uUS91iBwl1eO6omnjuae9ZjU4ew6cabenQ8vjbXWx7izdvOO2Chz1nIVipavihEsXR450iVB1NCjahSjsz1D/mbPUiPLWoPNcyIlXcbnmfU42E+DCSnbbrkv8nlVnuR5an9V1GDLVOYv/5SlyyNHuiQ10M3sBklvl3TE3V8d8fgrJX1R0nmSPurunyy9lCna0r7WlKzNKVlO0rwTZeJCNku/xujrDNrNs/zuapPafry1tXxdHTnSJalT/83sjZJ+IemmmEB/iaSXS9oi6Ymsgd72Xyxqq6jOQqm8tWbyricSN9U775T4qLH+w7q+nADihbQWUB0KTf1397vMbGPC40ckHTGzbiw+0mFxnZ/XXH5O7HjpvPJMukqb9RhV04+bpGUWvzzsJDSbTTI6O8uTaXGufqDfHlVDH3rOxyX9IqmGbmZbJW2VpA0bNrzu0UcfzVncyVbH6oyjhmtPL8gxyiVpgo10YnPAtlvu6/xCZUAdWrM4l7vvlLRT6jW51PneIWhiuNy4tae8sx7pHAOKC2KUS5KQ2ue6tjpjSCNHgC4IYhx6nDb94G8Z2rQ6Y1mrJA7UsXofELoso1xulnSBpDWSfiLpY5KmJcndP2tmvy5pQdLzJf1KvRExZ7v7z5Net45RLk20OVetDVccbfpRizZ8HkCdio5yuSLl8R9LWj9m2SrV9ina42jDiIBxF1oqW1UrSAJdFXSTS9unaHdV2kJLdTVxZVnyAJgkQQd6m9qc61R2+/aouC/EpIWWiojbnhCvwIAigg70Sexoq6MjuM6FlpK2hyswYLXghy22oc25KlX/ClCcOhdaStoehjoCqwUf6KGq41eAktS10FLarwlJ7VuECmgKgd5R4/wKUNWqCNi0yVQhX4EBeRHoHZX0K0Az01ONNUOUHbA0qwDZBd0pGrK4Gveg4zeUjuBJ7NgGxpVptcUqsB56MW2arQmgPq1ZbRHloUMQwCgCvcPoEAQwjDZ0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIFID3cxuMLMjZvZgzONmZp8xs4Nmdr+ZnVd+MQEAabLU0G+UdHHC45dIOrP/Z6uk64sXCwCQV2qgu/tdkn6W8JTLJN3kPd+SNGtmLyurgACAbMpoQ5+T9NjQ7UP9+wAANSoj0C3iPo98otlWM1sws4WlpaUS3hoAMFBGoB+SdPrQ7fWSDkc90d13uvu8u8+vXbu2hLcGAAyUEei7Jb23P9rlfElPuvvjJbwuACCHk9OeYGY3S7pA0hozOyTpY5KmJcndPyvpG5IulXRQ0lOS/rCqwgIA4qUGurtfkfK4S/pAaSUCAIyFmaIAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACESmQDezi83sgJkdNLPtEY+/3MzuMLP7zew/zGx9+UUFACRJDXQzm5J0naRLJJ0t6QozO3vkaZ+UdJO7v0bSJyRdU3ZBAQDJstTQXy/poLv/wN2flvQVSZeNPOdsSXf0/31nxOMAgIplCfQ5SY8N3T7Uv2/YtyW9s//v35X0a2b24tEXMrOtZrZgZgtLS0vjlBcAECNLoFvEfT5y+08kvcnM9kt6k6RFScdO+E/uO9193t3n165dm7uwAIB4J2d4ziFJpw/dXi/p8PAT3P2wpMslycyeJ+md7v5kWYUEAKTLUkO/V9KZZnaGmZ0i6d2Sdg8/wczWmNngta6WdEO5xQQApEkNdHc/JumDkvZIekjSre7+HTP7hJm9o/+0CyQdMLPvSXqppL+uqLwAgBjmPtocXo/5+XlfWFho5L0BoKvMbJ+7z0c9xkxRAAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBAEOgAE4uSmCwAAXbRr/6Ku3XNAh5dXtG52RldddJa2bJprtEwEOgDktGv/oq6+7QGtHD0uSVpcXtHVtz0gSY2GOk0uAJDTtXsOPBPmAytHj+vaPQcaKlEPgQ4AOR1eXsl1f10IdADIad3sTK7760KgA0BOV110lmamp1bdNzM9pasuOquhEvXQKQoAOQ06PhnlAgAB2LJprvEAH0WTCwAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgcgU6GZ2sZkdMLODZrY94vENZnanme03s/vN7NLyiwoASJI6scjMpiRdJ+l3JB2SdK+Z7Xb37w497c8l3eru15vZ2ZK+IWljBeUFOquN62cjLFlq6K+XdNDdf+DuT0v6iqTLRp7jkp7f//cLJB0ur4hA9w3Wz15cXpHr2fWzd+1fbLpoCEiWqf9zkh4bun1I0htGnvNxSf9mZn8k6TRJb416ITPbKmmrJG3YsCFvWYHOSlo/m1p6u3XpyipLDd0i7vOR21dIutHd10u6VNKXzeyE13b3ne4+7+7za9euzV9aoKPaun42knXtyipLoB+SdPrQ7fU6sUnlSkm3SpK7/6ekUyWtKaOAQAjaun42krX1l4niZAn0eyWdaWZnmNkpkt4taffIc34k6S2SZGavUi/Ql8os6CTbtX9Rm3fs1Rnbv67NO/a2tnaAeG1dPxvJunZlldqG7u7HzOyDkvZImpJ0g7t/x8w+IWnB3XdL+oikz5nZNvWaY/7A3UebZTCGtv4YLfJp6/rZSLZudkaLEeHd1israyp35+fnfWFhoZH37pLNO/ZGHlBzszO6e/uFDZQImByjFSqpd2V1zeXnNPZlbGb73H0+6jF+4KLlunbJB4Ska1dWBHrLde2SDwhNG3+ZKA5rubQcnWlAt9U5qIEaest17ZIPwLPqHtRAoHdAly75ADyr7hnCBDoAlCBqiYC6BzUQ6ABQUFzTyuxzp/XEU0dPeH5VgxoIdAAoKK5pZeXocZlWL35V5aAGAh0ACkpqQnHpmVCfq3hQA4EOAAXFzRcZGIR51bO7GYcOAAVFzRcZVcfsbmroAFDQ8HyRuJp6HbO7qaEDQAm2bJrT3dsv1KffdW5js7upoQNAiZqc3U2gA0DJmprdTZMLAASCQAeAQBDoABAIAh0AAkGgA0AgGvuRaDNbkvTomP99jaSfllicrpjE7Z7EbZYmc7sncZul/Nv9cndfG/VAY4FehJktxP3qdcgmcbsncZulydzuSdxmqdztpskFAAJBoANAILoa6DubLkBDJnG7J3Gbpcnc7kncZqnE7e5kGzoA4ERdraEDAEYQ6AAQiM4FupldbGYHzOygmW1vujxVMLPTzexOM3vIzL5jZh/q3/8iM/t3M/t+/+8XNl3WKpjZlJntN7Pb+7fPMLN7+tt9i5md0nQZy2Rms2b2VTN7uL/Pf3MS9rWZbesf3w+a2c1mdmqI+9rMbjCzI2b24NB9kfvXej7Tz7f7zey8PO/VqUA3sylJ10m6RNLZkq4ws7ObLVUljkn6iLu/StL5kj7Q387tku5w9zMl3dG/HaIPSXpo6PbfSPpUf7ufkHRlI6Wqzt9J+ld3f6Wk16q37UHvazObk/THkubd/dWSpiS9W2Hu6xslXTxyX9z+vUTSmf0/WyVdn+eNOhXokl4v6aC7/8Ddn5b0FUmXNVym0rn74+7+3/1//696J/icetv6pf7TviRpSzMlrI6ZrZf0Nkmf7982SRdK+mr/KUFtt5k9X9IbJX1Bktz9aXdf1gTsa/V+j2HGzE6W9FxJjyvAfe3ud0n62cjdcfv3Mkk3ec+3JM2a2cuyvlfXAn1O0mNDtw/17wuWmW2UtEnSPZJe6u6PS73Ql/SS5kpWmU9L+lNJv+rffrGkZXc/1r8d2j5/haQlSV/sNzN93sxOU+D72t0XJX1S0o/UC/InJe1T2Pt6WNz+LZRxXQt0i7gv2HGXZvY8SV+T9GF3/3nT5amamb1d0hF33zd8d8RTQ9rnJ0s6T9L17r5J0i8VWPNKlH6b8WWSzpC0TtJp6jU3jAppX2dR6HjvWqAfknT60O31kg43VJZKmdm0emH+D+5+W//unwwuv/p/H2mqfBXZLOkdZvZD9ZrTLlSvxj7bvyyXwtvnhyQdcvd7+re/ql7Ah76v3yrpEXdfcvejkm6T9FsKe18Pi9u/hTKua4F+r6Qz+z3hp6jXibK74TKVrt9u/AVJD7n73w49tFvS+/r/fp+kf667bFVy96vdfb27b1Rv3+519/dIulPS7/WfFtR2u/uPJT1mZoOfhH+LpO8q8H2tXlPL+Wb23P7xPtjuYPf1iLj9u1vSe/ujXc6X9OSgaSYTd+/UH0mXSvqepP+R9NGmy1PRNv62epdZ90u6r//nUvXak++Q9P3+3y9quqwVfgYXSLq9/+9XSPovSQcl/aOk5zRdvpK39VxJC/39vUvSCydhX0v6S0kPS3pQ0pclPSfEfS3pZvX6CY6qVwO/Mm7/qtfkcl0/3x5QbxRQ5vdi6j8ABKJrTS4AgBgEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAjE/wOEmW6CL6gdZgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot((params/bse)[idx_sort_aic[:k_best], 1], 'o')"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0xf829495308>]"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWc0lEQVR4nO3df4xdZZ3H8fdnhwFvWeUijEqHSiVLhlUrFm5o0Y0Bca1UAt3KxpIlCqvp4rL+yqbGxsRE46YxNausbGi6uCr+QFcotWvQ6ipGNxF2b1u0rDIrItrOoFzRKWqvy3T87h/3FG9v79w591fvzDOfVzKZc5/znHO+5z6nn9w5c6aPIgIzM1v4/mjQBZiZWW840M3MEuFANzNLhAPdzCwRDnQzs0ScNKgDn3nmmbF8+fJBHd7MbEHas2fPLyJipNm6gQX68uXLKZfLgzq8mdmCJOkns63zLRczs0Q40M3MEuFANzNLhAPdzCwRDnQzs0TkespF0juBNwMB7AduiIjf1a0/BbgduAh4Anh9RDza82rNemjnvgm27h5ncqrK0mKBTWvGWLdytGf9W21T335aYRgJpg5Pz7rc7rZ5+jeeQ7s19Xt5abHAZeePcO9DlRNe02zj28k10Mvt56K5/rdFSaPAfwIvjIiqpH8D7omIT9T1+VvgJRFxo6QNwF9ExOtb7bdUKoUfW7RB2blvgs079lOdnnm6rTA8xJb1K5r+A2u3f6ttXnfRKHftmTimPY92t83T/+g5AMfVutg1jm8n10C9brc/StKeiCg1W5f3lstJQEHSScASYLJh/dXAJ7PlO4HLJSl3hWYn2Nbd48eFV3V6hq27x3vSv9U2d9x/oKPgbHfbPP2PnkOzWhe7xvHt5Bqo1+32ecwZ6BExAXwI+CnwGHAoIr7a0G0UOJD1PwIcAs5o3JekjZLKksqVSqXb2s06NjlV7Wt7q3UzXcxB0O62efpPTlVbnsdiVv++dHIN5OnXy/d+zkCXdDq1T+AvAJYCp0q6rrFbk02Pu5IiYntElCKiNDLS9C9XzU6IpcVCX9tbrRvq4ofXdrfN039psdDyPBaz+velk2sgT79evvd5brm8CvhxRFQiYhrYAbysoc9BYBlAdlvmNOCXPavSrMc2rRmjMDx0TFtheIhNa8Z60r/VNteuWnZcex7tbpun/9FzaFbrYtc4vp1cA/W63T6PPE+5/BRYLWkJUAUuBxp/m7kLeCPwHeAa4Bvhue1sHjv6S6i8Txy023+ubUrnPLvjp1zybJunf7Nz8FMus783nVwD3V5D7ZrzKRcASe8DXg8cAfZRe4TxPUA5InZJegbwKWAltU/mGyLikVb79FMuZmbta/WUS65A7wcHuplZ+3rx2KKZmc1zDnQzs0Q40M3MEuFANzNLhAPdzCwRDnQzs0Q40M3MEuFANzNLhAPdzCwRDnQzs0Q40M3MEuFANzNLhAPdzCwRDnQzs0Q40M3MEuFANzNLRJ5JosckPVD39aSkdzT0uVTSobo+7+1fyWZm1sycc4pGxDjwUgBJQ8AEcHeTrt+OiCt7W56ZmeXV7i2Xy4EfRcRP+lGMmZl1rt1A3wDcMcu6SyR9V9KXJb2oWQdJGyWVJZUrlUqbhzYzs1ZyB7qkk4GrgC80Wb0XOCciLgA+Cuxsto+I2B4RpYgojYyMdFKvmZnNop1P6FcAeyPi540rIuLJiPhNtnwPMCzpzB7VaGZmObQT6Ncyy+0WSc+TpGz54my/T3RfnpmZ5TXnUy4AkpYAfw78TV3bjQARsQ24BniLpCNAFdgQEdH7cs3MbDa5Aj0iDgNnNLRtq1u+Bbilt6WZmVk7/JeiZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIuac4ELSGPD5uqZzgfdGxEfq+gi4GVgLHAauj4i9Pa513ti5b4Ktu8eZnKpyWmEYCaYOT8+6vLRY4LLzR7j3oUrubRq337RmjHUrR9s+9omo70Qu97LWdt/XPP0b66vfplfXXJ595u3f7n5tflM7M8VJGgImgFUR8ZO69rXAW6kF+irg5ohY1WpfpVIpyuVyR0UP0s59E2zesZ/q9MwJPW5heIjXXTTKXXsmTvixU9bu+9rJOBSGh9iyfkXHQdnsmmu1z7z9292vzQ+S9kREqdm6dm+5XA78qD7MM1cDt0fNfUBR0lkd1Drvbd09PpBArU7PcMf9BxzmPdbu+9rJOFSnZ9i6e7zTEptec632mbd/u/u1+a/dQN8A3NGkfRQ4UPf6YNZ2DEkbJZUllSuVSpuHnh8mp6oDO/aM593ui3bf107GoZvrZrZtu21vd3ub/3IHuqSTgauALzRb3aTtuKs+IrZHRCkiSiMjI/mrnEeWFgsDO/aQmr3N1q1239dOxqGb62a2bbttb3d7m//a+YR+BbA3In7eZN1BYFnd67OByW4Km682rRmjMDx0wo9bGB7i2lXLBnLslLX7vnYyDoXhITatGeu0xKbXXKt95u3f7n5t/msn0K+l+e0WgF3AG1SzGjgUEY91Xd08tG7lKFvWr2C0WEBAsTDM6UuGWy6PFgtct/r5bW3TuP2W9Sv4wLoVbR/7RNR3Ipd7WWu772ue/o31Hd2mm18yNl5zc+0zb/9292vzX66nXCQtoXaP/NyIOJS13QgQEduyxxZvAV5D7bHFGyKi5SMsC/UpFzOzQWr1lMucz6EDRMRh4IyGtm11ywHc1E2RZmbWHf+lqJlZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSUiV6BLKkq6U9JDkn4g6ZKG9ZdKOiTpgezrvf0p18zMZpNrggvgZuArEXFNNln0kiZ9vh0RV/auNDMza8ecgS7pWcArgOsBIuIp4Kn+lmVmZu3Kc8vlXKACfFzSPkm3STq1Sb9LJH1X0pclvajZjiRtlFSWVK5UKt3UbWZmDfIE+knAhcCtEbES+C3w7oY+e4FzIuIC4KPAzmY7iojtEVGKiNLIyEgXZZuZWaM8gX4QOBgR92ev76QW8E+LiCcj4jfZ8j3AsKQze1qpmZm1NGegR8TPgAOSxrKmy4Hv1/eR9DxJypYvzvb7RI9rNTOzFvI+5fJW4DPZEy6PADdIuhEgIrYB1wBvkXQEqAIbIiL6UbCZmTWnQeVuqVSKcrk8kGObmS1UkvZERKnZOv+lqJlZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlohcMxZJKgK3AS8GAvjriPhO3XoBNwNrgcPA9RGxt/fl/sHOfRNs3T3O5FSVpcUCl50/wr0PVZicqnJaYRgJpg5P92V5abHApjVjrFs52s9TNDNrS94p6G4GvhIR12TT0C1pWH8FcF72tQq4NfveFzv3TbB5x36q0zMATExV+fR9P316/VR1uq/LE1NVNu/YD+BQN7N5Y85bLpKeBbwC+BhARDwVEVMN3a4Gbo+a+4CipLN6Xm1m6+7xp8N8UKrTM2zdPT7QGszM6uW5h34uUAE+LmmfpNskndrQZxQ4UPf6YNZ2DEkbJZUllSuVSsdFT05VO962l+ZLHWZmkC/QTwIuBG6NiJXAb4F3N/RRk+2Om306IrZHRCkiSiMjI20Xe9TSYqHjbXtpvtRhZgb5Av0gcDAi7s9e30kt4Bv7LKt7fTYw2X15zW1aM0ZheKhfu8+lMDzEpjVjA63BzKzenIEeET8DDkg6ml6XA99v6LYLeINqVgOHIuKx3pb6B+tWjrJl/QpGiwUEjBYLXLf6+U+/LhaGOX3JcN+WR4sFtqxf4V+Imtm8kvcpl7cCn8mecHkEuEHSjQARsQ24h9ojiw9Te2zxhj7Ueox1K0cdqGZmdXIFekQ8AJQamrfVrQ/gph7WZWZmbfJfipqZJcKBbmaWCAe6mVkiHOhmZolwoJuZJcKBbmaWCAe6mVkiHOhmZolwoJuZJcKBbmaWCAe6mVkiHOhmZolwoJuZJcKBbmaWCAe6mVkiHOhmZonINcGFpEeBXwMzwJGIKDWsvxT4IvDjrGlHRLy/d2Wamdlc8k5BB3BZRPyixfpvR8SV3RZkZmad8S0XM7NE5A30AL4qaY+kjbP0uUTSdyV9WdKLmnWQtFFSWVK5Uql0VLCZmTWX95bLyyNiUtJzgK9JeigivlW3fi9wTkT8RtJaYCdwXuNOImI7sB2gVCpFl7WbmVmdXJ/QI2Iy+/44cDdwccP6JyPiN9nyPcCwpDN7XKuZmbUwZ6BLOlXSM48uA68GHmzo8zxJypYvzvb7RO/LNTOz2eS55fJc4O4sr08CPhsRX5F0I0BEbAOuAd4i6QhQBTZEhG+pmJmdQHMGekQ8AlzQpH1b3fItwC29Lc3MzNrhxxbNzBLhQDczS4QD3cwsEQ50M7NEONDNzBLhQDczS4QD3cwsEQ50M7NEONDNzBLhQDczS4QD3cwsEQ50M7NEONDNzBLhQDczS4QD3cwsEbnmFJX0KPBrYAY4EhGlhvUCbgbWAoeB6yNib29LhZ37Jti6e5zJqSpLiwU2rRlj3crRXh/GrC98/S4c9WN1WmEYCaYOTx8zbo3jedn5I9z7UOW4bWbbvh+UZ2KhLNBLEfGLWdavBd5KLdBXATdHxKpW+yyVSlEul3MXunPfBJt37Kc6PfN0W2F4iC3rV/gfhc17vn4XjmZjVa8wPMTrLhrlrj0Ts/Zppdtxl7Sn8UP1Ub265XI1cHvU3AcUJZ3Vo30DsHX3+HFvXnV6hq27x3t5GLO+8PW7cDQbq3rV6RnuuP9AR2F+dPt+jXveQA/gq5L2SNrYZP0ocKDu9cGs7RiSNkoqSypXKpW2Cp2cqrbVbjaf+PpdOPKMyUyXUyb3a9zzBvrLI+JC4ArgJkmvaFivJtscd8YRsT0iShFRGhkZaavQpcVCW+1m84mv34Ujz5gMqVnk9fYYncgV6BExmX1/HLgbuLihy0FgWd3rs4HJXhR41KY1YxSGh45pKwwPsWnNWC8PY9YXvn4XjmZjVa8wPMS1q5a17NNKP8d9zkCXdKqkZx5dBl4NPNjQbRfwBtWsBg5FxGO9LHTdylG2rF/BaLGAgNFiwb9QsgXD1+/C0ThWxcIwpy8ZPmbcPrBuxXHjed3q5zfdptn2A3vKRdK51D6VQ+0xx89GxD9IuhEgIrZljy3eAryG2mOLN0REy0dY2n3KxczMWj/lMudz6BHxCHBBk/ZtdcsB3NRNkWZm1h3/paiZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpaI3IEuaUjSPklfarLuekkVSQ9kX2/ubZlmZjaXOWcsqvN24AfAs2ZZ//mI+LvuSzIzs07k+oQu6WzgtcBt/S3HzMw6lfeWy0eAdwG/b9HndZK+J+lOScuadZC0UVJZUrlSqbRbq5mZtTBnoEu6Eng8Iva06PbvwPKIeAnwH8Anm3WKiO0RUYqI0sjISEcFm5lZc3k+ob8cuErSo8DngFdK+nR9h4h4IiL+L3v5L8BFPa3SzMzmNGegR8TmiDg7IpYDG4BvRMR19X0knVX38ipqvzw1M7MTqJ2nXI4h6f1AOSJ2AW+TdBVwBPglcH1vyjMzs7wUEQM5cKlUinK5PJBjm5ktVJL2RESp2Tr/paiZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpaI3DMWSRoCysBERFzZsO4U4HZqc4k+Abw+Ih7tYZ2L2s59E2zdPc7kVJWlxQKb1oyxbuXooMsys3mmnU/ob2f2uULfBPwqIv4E+DDwwW4Ls5qd+ybYvGM/E1NVApiYqrJ5x3527psYdGlmNs/kCnRJZwOvBW6bpcvVwCez5TuByyWp+/Js6+5xqtMzx7RVp2fYunt8QBWZ2XyV9xP6R4B3Ab+fZf0ocAAgIo4Ah4AzGjtJ2iipLKlcqVQ6KHfxmZyqttVuZovXnIEu6Urg8YjY06pbk7bjZp+OiO0RUYqI0sjISBtlLl5Li4W22s1s8crzCf3lwFWSHgU+B7xS0qcb+hwElgFIOgk4DfhlD+tctDatGaMwPHRMW2F4iE1rxgZUkZnNV3MGekRsjoizI2I5sAH4RkRc19BtF/DGbPmarM9xn9CtfetWjrJl/QpGiwUEjBYLbFm/wk+5mNlxcj+22EjS+4FyROwCPgZ8StLD1D6Zb+hRfUYt1B3gZjaXtgI9Ir4JfDNbfm9d+++Av+xlYWZm1h7/paiZWSIc6GZmiXCgm5klwoFuZpYIDerpQkkV4Ccdbn4m8IselrNQLMbzXoznDIvzvBfjOUP7531ORDT9y8yBBXo3JJUjojToOk60xXjei/GcYXGe92I8Z+jtefuWi5lZIhzoZmaJWKiBvn3QBQzIYjzvxXjOsDjPezGeM/TwvBfkPXQzMzveQv2EbmZmDRzoZmaJWHCBLuk1ksYlPSzp3YOupx8kLZN0r6QfSPofSW/P2p8t6WuSfph9P33QtfaDpCFJ+yR9KXv9Akn3Z+f9eUknD7rGXpJUlHSnpIeyMb9kMYy1pHdm1/eDku6Q9IwUx1rSv0p6XNKDdW1Nx1c1/5Tl2/ckXdjOsRZUoEsaAv4ZuAJ4IXCtpBcOtqq+OAL8fUT8KbAauCk7z3cDX4+I84CvZ69T1Dgh+QeBD2fn/Stqk5Kn5GbgKxFxPnABtXNPeqwljQJvA0oR8WJgiNp/u53iWH8CeE1D22zjewVwXva1Ebi1nQMtqEAHLgYejohHIuIpajMoXT3gmnouIh6LiL3Z8q+p/QMf5djJuD8JrBtMhf3TOCF5Ntn4K6lNPg6JnbekZwGvoDanABHxVERMsQjGmtp/313IZjlbAjxGgmMdEd/i+BncZhvfq4Hbo+Y+oCjprLzHWmiB/vRk1JmDWVuyJC0HVgL3A8+NiMegFvrAcwZXWd80Tkh+BjCVTT4O6Y35uUAF+Hh2m+k2SaeS+FhHxATwIeCn1IL8ELCHtMe63mzj21XGLbRAzzUZdSok/TFwF/COiHhy0PX02ywTkqc+5icBFwK3RsRK4Lckdnulmeye8dXAC4ClwKnUbjc0Smms8+jqel9ogf70ZNSZs4HJAdXSV5KGqYX5ZyJiR9b886M/fmXfHx9UfX1y3ITk1D6xF7MfyyG9MT8IHIyI+7PXd1IL+NTH+lXAjyOiEhHTwA7gZaQ91vVmG9+uMm6hBfp/A+dlvwk/mdovUXYNuKaey+4bfwz4QUT8Y92q+sm43wh88UTX1k+zTEj+V8C91CYfh8TOOyJ+BhyQNJY1XQ58n8THmtqtltWSlmTX+9HzTnasG8w2vruAN2RPu6wGDh29NZNLRCyoL2At8L/Aj4D3DLqePp3jn1H7Met7wAPZ11pq95O/Dvww+/7sQdfax/fgUuBL2fK5wH8BDwNfAE4ZdH09PteXAuVsvHcCpy+GsQbeBzwEPAh8CjglxbEG7qD2e4Jpap/A3zTb+FK75fLPWb7tp/YUUO5j+U//zcwSsdBuuZiZ2Swc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5kl4v8B4VTC00lqb6IAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot((nobs - df_resid)[idx_sort_aic[:k_best]], 'o')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"index in sorted array with small params value"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"20"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.argmin(params[idx_sort_aic[:80], 1])"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([4124, 4274, 705, 319, 792, 1453, 4733, 4275, 4715, 5558, 2230,\n",
" 1432, 689, 1298, 4714, 4264, 4131, 4286, 1647, 847, 1672],\n",
" dtype=int64)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.argsort(idx_sort_aic)[:20+1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"index of model with small params value in original model sequence"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"124"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"idx_sort_aic[20]"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"df_summ = res_all0.sorted_frame()"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"exog_idx [0, 1, 3, 8]\n",
"aic 293.801\n",
"Name: 20, dtype: object"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_summ.iloc[20]"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(293.80101229140456, array([1.03476655, 0.08190939]))"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res_all0.aic[124], params[124]"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[4, array([ True, True, False, True, False, False, False, False, True,\n",
" False, False, False, False, False, False]), 102.03370934823904, 293.80101229140456]"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res_all0.res_all[124]"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([290.93887686, 291.49160103, 292.39084523, 292.42139933,\n",
" 292.49144156, 292.61745297, 292.67275953, 292.75913778,\n",
" 292.89747028, 292.90781 , 292.93855279, 293.01250351,\n",
" 293.02511464, 293.16586405, 293.23282997, 293.2569139 ,\n",
" 293.27831641, 293.32635662, 293.445069 , 293.48736274,\n",
" 293.80101229, 293.99163866, 294.04939289, 294.11036986,\n",
" 294.11841962, 294.17896341, 294.22853142, 294.23066083,\n",
" 294.23776358, 294.23946088, 294.25376474, 294.26523045,\n",
" 294.32214283, 294.32675852, 294.35873889, 294.37538128,\n",
" 294.38068857, 294.39444662, 294.3986518 , 294.40839955,\n",
" 294.42055476, 294.45667116, 294.47376241, 294.48703023,\n",
" 294.48761783, 294.58700314, 294.60835245, 294.614979 ,\n",
" 294.61592131, 294.62090463])"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res_all0.aic[idx_sort_aic[:50]]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Add parameter of interest to sorted results dataframe"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>exog_idx</th>\n",
" <th>aic</th>\n",
" <th>par1</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>[0, 1, 3, 8, 12]</td>\n",
" <td>290.938877</td>\n",
" <td>0.117691</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>[0, 1, 3, 6, 8, 12]</td>\n",
" <td>291.491601</td>\n",
" <td>0.120413</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>[0, 1, 3, 7, 8, 12]</td>\n",
" <td>292.390845</td>\n",
" <td>0.115739</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>[0, 1, 2, 3, 8, 12]</td>\n",
" <td>292.421399</td>\n",
" <td>0.115959</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>[0, 1, 3, 5, 8, 12]</td>\n",
" <td>292.491442</td>\n",
" <td>0.113315</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>[0, 1, 3, 8, 11, 12]</td>\n",
" <td>292.617453</td>\n",
" <td>0.121442</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>[0, 1, 3, 8, 12, 13]</td>\n",
" <td>292.672760</td>\n",
" <td>0.117442</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>[0, 1, 3, 8, 10, 12]</td>\n",
" <td>292.759138</td>\n",
" <td>0.117169</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>[0, 1, 3, 8, 12, 14]</td>\n",
" <td>292.897470</td>\n",
" <td>0.117100</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>[0, 1, 3, 8, 9, 12]</td>\n",
" <td>292.907810</td>\n",
" <td>0.119063</td>\n",
" </tr>\n",
" <tr>\n",
" <th>10</th>\n",
" <td>[0, 1, 3, 4, 8, 12]</td>\n",
" <td>292.938553</td>\n",
" <td>0.117864</td>\n",
" </tr>\n",
" <tr>\n",
" <th>11</th>\n",
" <td>[0, 1, 3, 6, 8, 11, 12]</td>\n",
" <td>293.012504</td>\n",
" <td>0.125133</td>\n",
" </tr>\n",
" <tr>\n",
" <th>12</th>\n",
" <td>[0, 1, 3, 6, 7, 8, 12]</td>\n",
" <td>293.025115</td>\n",
" <td>0.118542</td>\n",
" </tr>\n",
" <tr>\n",
" <th>13</th>\n",
" <td>[0, 1, 3, 6, 8, 12, 14]</td>\n",
" <td>293.165864</td>\n",
" <td>0.119080</td>\n",
" </tr>\n",
" <tr>\n",
" <th>14</th>\n",
" <td>[0, 1, 3, 6, 8, 12, 13]</td>\n",
" <td>293.232830</td>\n",
" <td>0.120160</td>\n",
" </tr>\n",
" <tr>\n",
" <th>15</th>\n",
" <td>[0, 1, 2, 3, 6, 8, 12]</td>\n",
" <td>293.256914</td>\n",
" <td>0.118998</td>\n",
" </tr>\n",
" <tr>\n",
" <th>16</th>\n",
" <td>[0, 1, 3, 6, 8, 10, 12]</td>\n",
" <td>293.278316</td>\n",
" <td>0.119879</td>\n",
" </tr>\n",
" <tr>\n",
" <th>17</th>\n",
" <td>[0, 1, 3, 5, 6, 8, 12]</td>\n",
" <td>293.326357</td>\n",
" <td>0.117482</td>\n",
" </tr>\n",
" <tr>\n",
" <th>18</th>\n",
" <td>[0, 1, 3, 6, 8, 9, 12]</td>\n",
" <td>293.445069</td>\n",
" <td>0.122097</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19</th>\n",
" <td>[0, 1, 3, 4, 6, 8, 12]</td>\n",
" <td>293.487363</td>\n",
" <td>0.121043</td>\n",
" </tr>\n",
" <tr>\n",
" <th>20</th>\n",
" <td>[0, 1, 3, 8]</td>\n",
" <td>293.801012</td>\n",
" <td>0.081909</td>\n",
" </tr>\n",
" <tr>\n",
" <th>21</th>\n",
" <td>[0, 1, 2, 3, 7, 8, 12]</td>\n",
" <td>293.991639</td>\n",
" <td>0.114422</td>\n",
" </tr>\n",
" <tr>\n",
" <th>22</th>\n",
" <td>[0, 1, 2, 3, 5, 8, 12]</td>\n",
" <td>294.049393</td>\n",
" <td>0.112091</td>\n",
" </tr>\n",
" <tr>\n",
" <th>23</th>\n",
" <td>[0, 1, 2, 3, 8, 12, 13]</td>\n",
" <td>294.110370</td>\n",
" <td>0.115617</td>\n",
" </tr>\n",
" <tr>\n",
" <th>24</th>\n",
" <td>[0, 1, 3, 5, 7, 8, 12]</td>\n",
" <td>294.118420</td>\n",
" <td>0.112555</td>\n",
" </tr>\n",
" <tr>\n",
" <th>25</th>\n",
" <td>[0, 1, 3, 7, 8, 11, 12]</td>\n",
" <td>294.178963</td>\n",
" <td>0.119002</td>\n",
" </tr>\n",
" <tr>\n",
" <th>26</th>\n",
" <td>[0, 1, 3, 5, 8, 11, 12]</td>\n",
" <td>294.228531</td>\n",
" <td>0.117000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>27</th>\n",
" <td>[0, 1, 2, 3, 8, 11, 12]</td>\n",
" <td>294.230661</td>\n",
" <td>0.119106</td>\n",
" </tr>\n",
" <tr>\n",
" <th>28</th>\n",
" <td>[0, 1, 3, 8, 11, 12, 13]</td>\n",
" <td>294.237764</td>\n",
" <td>0.121806</td>\n",
" </tr>\n",
" <tr>\n",
" <th>29</th>\n",
" <td>[0, 1, 3, 7, 8, 12, 13]</td>\n",
" <td>294.239461</td>\n",
" <td>0.115738</td>\n",
" </tr>\n",
" <tr>\n",
" <th>30</th>\n",
" <td>[0, 1, 2, 3, 8, 10, 12]</td>\n",
" <td>294.253765</td>\n",
" <td>0.115477</td>\n",
" </tr>\n",
" <tr>\n",
" <th>31</th>\n",
" <td>[0, 1, 3, 7, 8, 10, 12]</td>\n",
" <td>294.265230</td>\n",
" <td>0.115394</td>\n",
" </tr>\n",
" <tr>\n",
" <th>32</th>\n",
" <td>[0, 1, 3, 5, 8, 12, 13]</td>\n",
" <td>294.322143</td>\n",
" <td>0.113566</td>\n",
" </tr>\n",
" <tr>\n",
" <th>33</th>\n",
" <td>[0, 1, 3, 7, 8, 9, 12]</td>\n",
" <td>294.326759</td>\n",
" <td>0.117647</td>\n",
" </tr>\n",
" <tr>\n",
" <th>34</th>\n",
" <td>[0, 1, 3, 5, 8, 10, 12]</td>\n",
" <td>294.358739</td>\n",
" <td>0.113087</td>\n",
" </tr>\n",
" <tr>\n",
" <th>35</th>\n",
" <td>[0, 1, 3, 7, 8, 12, 14]</td>\n",
" <td>294.375381</td>\n",
" <td>0.115413</td>\n",
" </tr>\n",
" <tr>\n",
" <th>36</th>\n",
" <td>[0, 1, 3, 4, 7, 8, 12]</td>\n",
" <td>294.380689</td>\n",
" <td>0.116686</td>\n",
" </tr>\n",
" <tr>\n",
" <th>37</th>\n",
" <td>[0, 1, 2, 3, 8, 12, 14]</td>\n",
" <td>294.394447</td>\n",
" <td>0.115505</td>\n",
" </tr>\n",
" <tr>\n",
" <th>38</th>\n",
" <td>[0, 1, 2, 3, 8, 9, 12]</td>\n",
" <td>294.398652</td>\n",
" <td>0.117144</td>\n",
" </tr>\n",
" <tr>\n",
" <th>39</th>\n",
" <td>[0, 1, 3, 5, 8, 9, 12]</td>\n",
" <td>294.408400</td>\n",
" <td>0.115276</td>\n",
" </tr>\n",
" <tr>\n",
" <th>40</th>\n",
" <td>[0, 1, 2, 3, 4, 8, 12]</td>\n",
" <td>294.420555</td>\n",
" <td>0.116237</td>\n",
" </tr>\n",
" <tr>\n",
" <th>41</th>\n",
" <td>[0, 1, 3, 5, 8, 12, 14]</td>\n",
" <td>294.456671</td>\n",
" <td>0.112806</td>\n",
" </tr>\n",
" <tr>\n",
" <th>42</th>\n",
" <td>[0, 1, 3, 8, 10, 11, 12]</td>\n",
" <td>294.473762</td>\n",
" <td>0.120768</td>\n",
" </tr>\n",
" <tr>\n",
" <th>43</th>\n",
" <td>[0, 1, 3, 4, 5, 8, 12]</td>\n",
" <td>294.487030</td>\n",
" <td>0.113925</td>\n",
" </tr>\n",
" <tr>\n",
" <th>44</th>\n",
" <td>[0, 1, 3, 8, 10, 12, 13]</td>\n",
" <td>294.487618</td>\n",
" <td>0.116911</td>\n",
" </tr>\n",
" <tr>\n",
" <th>45</th>\n",
" <td>[0, 1, 3, 8, 9, 11, 12]</td>\n",
" <td>294.587003</td>\n",
" <td>0.122795</td>\n",
" </tr>\n",
" <tr>\n",
" <th>46</th>\n",
" <td>[0, 1, 3, 8, 11, 12, 14]</td>\n",
" <td>294.608352</td>\n",
" <td>0.121036</td>\n",
" </tr>\n",
" <tr>\n",
" <th>47</th>\n",
" <td>[0, 1, 3, 4, 8, 11, 12]</td>\n",
" <td>294.614979</td>\n",
" <td>0.121941</td>\n",
" </tr>\n",
" <tr>\n",
" <th>48</th>\n",
" <td>[0, 1, 3, 6, 8, 11, 12, 13]</td>\n",
" <td>294.615921</td>\n",
" <td>0.125514</td>\n",
" </tr>\n",
" <tr>\n",
" <th>49</th>\n",
" <td>[0, 1, 3, 8, 12, 13, 14]</td>\n",
" <td>294.620905</td>\n",
" <td>0.116776</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" exog_idx aic par1\n",
"0 [0, 1, 3, 8, 12] 290.938877 0.117691\n",
"1 [0, 1, 3, 6, 8, 12] 291.491601 0.120413\n",
"2 [0, 1, 3, 7, 8, 12] 292.390845 0.115739\n",
"3 [0, 1, 2, 3, 8, 12] 292.421399 0.115959\n",
"4 [0, 1, 3, 5, 8, 12] 292.491442 0.113315\n",
"5 [0, 1, 3, 8, 11, 12] 292.617453 0.121442\n",
"6 [0, 1, 3, 8, 12, 13] 292.672760 0.117442\n",
"7 [0, 1, 3, 8, 10, 12] 292.759138 0.117169\n",
"8 [0, 1, 3, 8, 12, 14] 292.897470 0.117100\n",
"9 [0, 1, 3, 8, 9, 12] 292.907810 0.119063\n",
"10 [0, 1, 3, 4, 8, 12] 292.938553 0.117864\n",
"11 [0, 1, 3, 6, 8, 11, 12] 293.012504 0.125133\n",
"12 [0, 1, 3, 6, 7, 8, 12] 293.025115 0.118542\n",
"13 [0, 1, 3, 6, 8, 12, 14] 293.165864 0.119080\n",
"14 [0, 1, 3, 6, 8, 12, 13] 293.232830 0.120160\n",
"15 [0, 1, 2, 3, 6, 8, 12] 293.256914 0.118998\n",
"16 [0, 1, 3, 6, 8, 10, 12] 293.278316 0.119879\n",
"17 [0, 1, 3, 5, 6, 8, 12] 293.326357 0.117482\n",
"18 [0, 1, 3, 6, 8, 9, 12] 293.445069 0.122097\n",
"19 [0, 1, 3, 4, 6, 8, 12] 293.487363 0.121043\n",
"20 [0, 1, 3, 8] 293.801012 0.081909\n",
"21 [0, 1, 2, 3, 7, 8, 12] 293.991639 0.114422\n",
"22 [0, 1, 2, 3, 5, 8, 12] 294.049393 0.112091\n",
"23 [0, 1, 2, 3, 8, 12, 13] 294.110370 0.115617\n",
"24 [0, 1, 3, 5, 7, 8, 12] 294.118420 0.112555\n",
"25 [0, 1, 3, 7, 8, 11, 12] 294.178963 0.119002\n",
"26 [0, 1, 3, 5, 8, 11, 12] 294.228531 0.117000\n",
"27 [0, 1, 2, 3, 8, 11, 12] 294.230661 0.119106\n",
"28 [0, 1, 3, 8, 11, 12, 13] 294.237764 0.121806\n",
"29 [0, 1, 3, 7, 8, 12, 13] 294.239461 0.115738\n",
"30 [0, 1, 2, 3, 8, 10, 12] 294.253765 0.115477\n",
"31 [0, 1, 3, 7, 8, 10, 12] 294.265230 0.115394\n",
"32 [0, 1, 3, 5, 8, 12, 13] 294.322143 0.113566\n",
"33 [0, 1, 3, 7, 8, 9, 12] 294.326759 0.117647\n",
"34 [0, 1, 3, 5, 8, 10, 12] 294.358739 0.113087\n",
"35 [0, 1, 3, 7, 8, 12, 14] 294.375381 0.115413\n",
"36 [0, 1, 3, 4, 7, 8, 12] 294.380689 0.116686\n",
"37 [0, 1, 2, 3, 8, 12, 14] 294.394447 0.115505\n",
"38 [0, 1, 2, 3, 8, 9, 12] 294.398652 0.117144\n",
"39 [0, 1, 3, 5, 8, 9, 12] 294.408400 0.115276\n",
"40 [0, 1, 2, 3, 4, 8, 12] 294.420555 0.116237\n",
"41 [0, 1, 3, 5, 8, 12, 14] 294.456671 0.112806\n",
"42 [0, 1, 3, 8, 10, 11, 12] 294.473762 0.120768\n",
"43 [0, 1, 3, 4, 5, 8, 12] 294.487030 0.113925\n",
"44 [0, 1, 3, 8, 10, 12, 13] 294.487618 0.116911\n",
"45 [0, 1, 3, 8, 9, 11, 12] 294.587003 0.122795\n",
"46 [0, 1, 3, 8, 11, 12, 14] 294.608352 0.121036\n",
"47 [0, 1, 3, 4, 8, 11, 12] 294.614979 0.121941\n",
"48 [0, 1, 3, 6, 8, 11, 12, 13] 294.615921 0.125514\n",
"49 [0, 1, 3, 8, 12, 13, 14] 294.620905 0.116776"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_summ['par1'] = params[idx_sort_aic, 1]\n",
"df_summ.iloc[:50]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"list models with small estimate of parameter of interest, sorted by aic\n",
"\n",
"Note, all models include explanatory variable 8. This indicates that variable 8 captures some of the positive effect of the variable of interest."
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>exog_idx</th>\n",
" <th>aic</th>\n",
" <th>par1</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>20</th>\n",
" <td>[0, 1, 3, 8]</td>\n",
" <td>293.801012</td>\n",
" <td>0.081909</td>\n",
" </tr>\n",
" <tr>\n",
" <th>57</th>\n",
" <td>[0, 1, 3, 6, 8]</td>\n",
" <td>294.791141</td>\n",
" <td>0.082995</td>\n",
" </tr>\n",
" <tr>\n",
" <th>93</th>\n",
" <td>[0, 1, 3, 8, 13]</td>\n",
" <td>295.285584</td>\n",
" <td>0.082727</td>\n",
" </tr>\n",
" <tr>\n",
" <th>95</th>\n",
" <td>[0, 1, 3, 8, 11]</td>\n",
" <td>295.334311</td>\n",
" <td>0.087203</td>\n",
" </tr>\n",
" <tr>\n",
" <th>96</th>\n",
" <td>[0, 1, 2, 3, 8]</td>\n",
" <td>295.342805</td>\n",
" <td>0.080112</td>\n",
" </tr>\n",
" <tr>\n",
" <th>97</th>\n",
" <td>[0, 1, 3, 7, 8]</td>\n",
" <td>295.381120</td>\n",
" <td>0.079800</td>\n",
" </tr>\n",
" <tr>\n",
" <th>98</th>\n",
" <td>[0, 1, 3, 8, 9]</td>\n",
" <td>295.398978</td>\n",
" <td>0.080992</td>\n",
" </tr>\n",
" <tr>\n",
" <th>100</th>\n",
" <td>[0, 1, 3, 5, 8]</td>\n",
" <td>295.464322</td>\n",
" <td>0.077706</td>\n",
" </tr>\n",
" <tr>\n",
" <th>102</th>\n",
" <td>[0, 1, 3, 8, 10]</td>\n",
" <td>295.567308</td>\n",
" <td>0.081554</td>\n",
" </tr>\n",
" <tr>\n",
" <th>103</th>\n",
" <td>[0, 1, 3, 8, 14]</td>\n",
" <td>295.595339</td>\n",
" <td>0.081446</td>\n",
" </tr>\n",
" <tr>\n",
" <th>106</th>\n",
" <td>[0, 1, 3, 4, 8]</td>\n",
" <td>295.800614</td>\n",
" <td>0.082107</td>\n",
" </tr>\n",
" <tr>\n",
" <th>140</th>\n",
" <td>[0, 1, 3, 6, 8, 11]</td>\n",
" <td>296.161910</td>\n",
" <td>0.089226</td>\n",
" </tr>\n",
" <tr>\n",
" <th>144</th>\n",
" <td>[0, 1, 3, 6, 8, 14]</td>\n",
" <td>296.180188</td>\n",
" <td>0.082420</td>\n",
" </tr>\n",
" <tr>\n",
" <th>163</th>\n",
" <td>[0, 1, 3, 6, 8, 13]</td>\n",
" <td>296.272639</td>\n",
" <td>0.083809</td>\n",
" </tr>\n",
" <tr>\n",
" <th>174</th>\n",
" <td>[0, 1, 3, 6, 8, 9]</td>\n",
" <td>296.387351</td>\n",
" <td>0.082079</td>\n",
" </tr>\n",
" <tr>\n",
" <th>181</th>\n",
" <td>[0, 1, 3, 6, 7, 8]</td>\n",
" <td>296.435048</td>\n",
" <td>0.081023</td>\n",
" </tr>\n",
" <tr>\n",
" <th>192</th>\n",
" <td>[0, 1, 3, 6, 8, 10]</td>\n",
" <td>296.523491</td>\n",
" <td>0.082634</td>\n",
" </tr>\n",
" <tr>\n",
" <th>195</th>\n",
" <td>[0, 1, 2, 3, 6, 8]</td>\n",
" <td>296.559356</td>\n",
" <td>0.081584</td>\n",
" </tr>\n",
" <tr>\n",
" <th>213</th>\n",
" <td>[0, 1, 3, 8, 11, 13]</td>\n",
" <td>296.636698</td>\n",
" <td>0.089167</td>\n",
" </tr>\n",
" <tr>\n",
" <th>216</th>\n",
" <td>[0, 1, 3, 5, 6, 8]</td>\n",
" <td>296.662078</td>\n",
" <td>0.080240</td>\n",
" </tr>\n",
" <tr>\n",
" <th>234</th>\n",
" <td>[0, 1, 2, 3, 8, 13]</td>\n",
" <td>296.767123</td>\n",
" <td>0.080863</td>\n",
" </tr>\n",
" <tr>\n",
" <th>237</th>\n",
" <td>[0, 1, 3, 4, 6, 8]</td>\n",
" <td>296.787634</td>\n",
" <td>0.083582</td>\n",
" </tr>\n",
" <tr>\n",
" <th>264</th>\n",
" <td>[0, 1, 2, 3, 8, 9]</td>\n",
" <td>296.904996</td>\n",
" <td>0.079090</td>\n",
" </tr>\n",
" <tr>\n",
" <th>273</th>\n",
" <td>[0, 1, 3, 8, 9, 13]</td>\n",
" <td>296.944944</td>\n",
" <td>0.081833</td>\n",
" </tr>\n",
" <tr>\n",
" <th>274</th>\n",
" <td>[0, 1, 3, 8, 9, 11]</td>\n",
" <td>296.946204</td>\n",
" <td>0.086215</td>\n",
" </tr>\n",
" <tr>\n",
" <th>282</th>\n",
" <td>[0, 1, 3, 7, 8, 13]</td>\n",
" <td>297.004175</td>\n",
" <td>0.080867</td>\n",
" </tr>\n",
" <tr>\n",
" <th>283</th>\n",
" <td>[0, 1, 2, 3, 7, 8]</td>\n",
" <td>297.020820</td>\n",
" <td>0.078446</td>\n",
" </tr>\n",
" <tr>\n",
" <th>285</th>\n",
" <td>[0, 1, 2, 3, 8, 11]</td>\n",
" <td>297.024460</td>\n",
" <td>0.084853</td>\n",
" </tr>\n",
" <tr>\n",
" <th>286</th>\n",
" <td>[0, 1, 3, 7, 8, 11]</td>\n",
" <td>297.028994</td>\n",
" <td>0.084736</td>\n",
" </tr>\n",
" <tr>\n",
" <th>290</th>\n",
" <td>[0, 1, 3, 7, 8, 9]</td>\n",
" <td>297.039983</td>\n",
" <td>0.079110</td>\n",
" </tr>\n",
" <tr>\n",
" <th>292</th>\n",
" <td>[0, 1, 3, 8, 10, 13]</td>\n",
" <td>297.045144</td>\n",
" <td>0.082372</td>\n",
" </tr>\n",
" <tr>\n",
" <th>294</th>\n",
" <td>[0, 1, 3, 8, 13, 14]</td>\n",
" <td>297.055920</td>\n",
" <td>0.082256</td>\n",
" </tr>\n",
" <tr>\n",
" <th>295</th>\n",
" <td>[0, 1, 3, 5, 8, 11]</td>\n",
" <td>297.057259</td>\n",
" <td>0.083045</td>\n",
" </tr>\n",
" <tr>\n",
" <th>296</th>\n",
" <td>[0, 1, 3, 5, 8, 13]</td>\n",
" <td>297.063289</td>\n",
" <td>0.079186</td>\n",
" </tr>\n",
" <tr>\n",
" <th>297</th>\n",
" <td>[0, 1, 2, 3, 5, 8]</td>\n",
" <td>297.068120</td>\n",
" <td>0.076431</td>\n",
" </tr>\n",
" <tr>\n",
" <th>304</th>\n",
" <td>[0, 1, 2, 3, 8, 10]</td>\n",
" <td>297.121803</td>\n",
" <td>0.079794</td>\n",
" </tr>\n",
" <tr>\n",
" <th>306</th>\n",
" <td>[0, 1, 3, 5, 8, 9]</td>\n",
" <td>297.150380</td>\n",
" <td>0.077453</td>\n",
" </tr>\n",
" <tr>\n",
" <th>307</th>\n",
" <td>[0, 1, 3, 8, 10, 11]</td>\n",
" <td>297.150796</td>\n",
" <td>0.086608</td>\n",
" </tr>\n",
" <tr>\n",
" <th>309</th>\n",
" <td>[0, 1, 3, 8, 9, 10]</td>\n",
" <td>297.156185</td>\n",
" <td>0.080622</td>\n",
" </tr>\n",
" <tr>\n",
" <th>311</th>\n",
" <td>[0, 1, 2, 3, 8, 14]</td>\n",
" <td>297.167600</td>\n",
" <td>0.079744</td>\n",
" </tr>\n",
" <tr>\n",
" <th>312</th>\n",
" <td>[0, 1, 3, 5, 7, 8]</td>\n",
" <td>297.177963</td>\n",
" <td>0.076787</td>\n",
" </tr>\n",
" <tr>\n",
" <th>314</th>\n",
" <td>[0, 1, 3, 7, 8, 10]</td>\n",
" <td>297.201493</td>\n",
" <td>0.079621</td>\n",
" </tr>\n",
" <tr>\n",
" <th>315</th>\n",
" <td>[0, 1, 3, 8, 9, 14]</td>\n",
" <td>297.206120</td>\n",
" <td>0.080560</td>\n",
" </tr>\n",
" <tr>\n",
" <th>319</th>\n",
" <td>[0, 1, 3]</td>\n",
" <td>297.222624</td>\n",
" <td>0.083677</td>\n",
" </tr>\n",
" <tr>\n",
" <th>320</th>\n",
" <td>[0, 1, 3, 8, 11, 14]</td>\n",
" <td>297.229051</td>\n",
" <td>0.086354</td>\n",
" </tr>\n",
" <tr>\n",
" <th>321</th>\n",
" <td>[0, 1, 3, 7, 8, 14]</td>\n",
" <td>297.229623</td>\n",
" <td>0.079532</td>\n",
" </tr>\n",
" <tr>\n",
" <th>323</th>\n",
" <td>[0, 1, 3, 4, 8, 13]</td>\n",
" <td>297.268742</td>\n",
" <td>0.084047</td>\n",
" </tr>\n",
" <tr>\n",
" <th>324</th>\n",
" <td>[0, 1, 3, 5, 8, 14]</td>\n",
" <td>297.269373</td>\n",
" <td>0.077327</td>\n",
" </tr>\n",
" <tr>\n",
" <th>325</th>\n",
" <td>[0, 1, 3, 5, 8, 10]</td>\n",
" <td>297.277023</td>\n",
" <td>0.077676</td>\n",
" </tr>\n",
" <tr>\n",
" <th>327</th>\n",
" <td>[0, 1, 3, 8, 10, 14]</td>\n",
" <td>297.329807</td>\n",
" <td>0.081032</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" exog_idx aic par1\n",
"20 [0, 1, 3, 8] 293.801012 0.081909\n",
"57 [0, 1, 3, 6, 8] 294.791141 0.082995\n",
"93 [0, 1, 3, 8, 13] 295.285584 0.082727\n",
"95 [0, 1, 3, 8, 11] 295.334311 0.087203\n",
"96 [0, 1, 2, 3, 8] 295.342805 0.080112\n",
"97 [0, 1, 3, 7, 8] 295.381120 0.079800\n",
"98 [0, 1, 3, 8, 9] 295.398978 0.080992\n",
"100 [0, 1, 3, 5, 8] 295.464322 0.077706\n",
"102 [0, 1, 3, 8, 10] 295.567308 0.081554\n",
"103 [0, 1, 3, 8, 14] 295.595339 0.081446\n",
"106 [0, 1, 3, 4, 8] 295.800614 0.082107\n",
"140 [0, 1, 3, 6, 8, 11] 296.161910 0.089226\n",
"144 [0, 1, 3, 6, 8, 14] 296.180188 0.082420\n",
"163 [0, 1, 3, 6, 8, 13] 296.272639 0.083809\n",
"174 [0, 1, 3, 6, 8, 9] 296.387351 0.082079\n",
"181 [0, 1, 3, 6, 7, 8] 296.435048 0.081023\n",
"192 [0, 1, 3, 6, 8, 10] 296.523491 0.082634\n",
"195 [0, 1, 2, 3, 6, 8] 296.559356 0.081584\n",
"213 [0, 1, 3, 8, 11, 13] 296.636698 0.089167\n",
"216 [0, 1, 3, 5, 6, 8] 296.662078 0.080240\n",
"234 [0, 1, 2, 3, 8, 13] 296.767123 0.080863\n",
"237 [0, 1, 3, 4, 6, 8] 296.787634 0.083582\n",
"264 [0, 1, 2, 3, 8, 9] 296.904996 0.079090\n",
"273 [0, 1, 3, 8, 9, 13] 296.944944 0.081833\n",
"274 [0, 1, 3, 8, 9, 11] 296.946204 0.086215\n",
"282 [0, 1, 3, 7, 8, 13] 297.004175 0.080867\n",
"283 [0, 1, 2, 3, 7, 8] 297.020820 0.078446\n",
"285 [0, 1, 2, 3, 8, 11] 297.024460 0.084853\n",
"286 [0, 1, 3, 7, 8, 11] 297.028994 0.084736\n",
"290 [0, 1, 3, 7, 8, 9] 297.039983 0.079110\n",
"292 [0, 1, 3, 8, 10, 13] 297.045144 0.082372\n",
"294 [0, 1, 3, 8, 13, 14] 297.055920 0.082256\n",
"295 [0, 1, 3, 5, 8, 11] 297.057259 0.083045\n",
"296 [0, 1, 3, 5, 8, 13] 297.063289 0.079186\n",
"297 [0, 1, 2, 3, 5, 8] 297.068120 0.076431\n",
"304 [0, 1, 2, 3, 8, 10] 297.121803 0.079794\n",
"306 [0, 1, 3, 5, 8, 9] 297.150380 0.077453\n",
"307 [0, 1, 3, 8, 10, 11] 297.150796 0.086608\n",
"309 [0, 1, 3, 8, 9, 10] 297.156185 0.080622\n",
"311 [0, 1, 2, 3, 8, 14] 297.167600 0.079744\n",
"312 [0, 1, 3, 5, 7, 8] 297.177963 0.076787\n",
"314 [0, 1, 3, 7, 8, 10] 297.201493 0.079621\n",
"315 [0, 1, 3, 8, 9, 14] 297.206120 0.080560\n",
"319 [0, 1, 3] 297.222624 0.083677\n",
"320 [0, 1, 3, 8, 11, 14] 297.229051 0.086354\n",
"321 [0, 1, 3, 7, 8, 14] 297.229623 0.079532\n",
"323 [0, 1, 3, 4, 8, 13] 297.268742 0.084047\n",
"324 [0, 1, 3, 5, 8, 14] 297.269373 0.077327\n",
"325 [0, 1, 3, 5, 8, 10] 297.277023 0.077676\n",
"327 [0, 1, 3, 8, 10, 14] 297.329807 0.081032"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_summ[df_summ[\"par1\"] < 0.09].iloc[:50]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Histogram of parameter estimates across model specification\n",
"\n",
"The histogram of the parameter of interest estimated with all subset models shows a tri-modal distribution.\n",
"If we only consider the best 100 models ranked by aic, then we observe a bimodal distribution without the upper part of all sets."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([ 26., 35., 99., 112., 143., 252., 252., 191., 189., 191., 181.,\n",
" 167., 112., 125., 181., 307., 354., 400., 497., 415., 363., 215.,\n",
" 185., 219., 199., 249., 222., 173., 69., 17., 4., 0., 0.,\n",
" 6., 21., 44., 108., 168., 246., 291., 301., 231., 218., 161.,\n",
" 116., 57., 43., 20., 11., 6.]),\n",
" array([0.05440995, 0.056015 , 0.05762006, 0.05922511, 0.06083017,\n",
" 0.06243522, 0.06404028, 0.06564533, 0.06725038, 0.06885544,\n",
" 0.07046049, 0.07206555, 0.0736706 , 0.07527566, 0.07688071,\n",
" 0.07848577, 0.08009082, 0.08169587, 0.08330093, 0.08490598,\n",
" 0.08651104, 0.08811609, 0.08972115, 0.0913262 , 0.09293126,\n",
" 0.09453631, 0.09614136, 0.09774642, 0.09935147, 0.10095653,\n",
" 0.10256158, 0.10416664, 0.10577169, 0.10737675, 0.1089818 ,\n",
" 0.11058685, 0.11219191, 0.11379696, 0.11540202, 0.11700707,\n",
" 0.11861213, 0.12021718, 0.12182224, 0.12342729, 0.12503235,\n",
" 0.1266374 , 0.12824245, 0.12984751, 0.13145256, 0.13305762,\n",
" 0.13466267]),\n",
" <a list of 50 Patch objects>)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAQgElEQVR4nO3dbYylZX3H8e9PVkDxYRdYCO6uDsY10bSKuiIJfbBQFcG6vICKWt1Skk1Tmtja1GLbxKbpC2iaUk0bmk3RLqYVkNayQSLFRat9wLqLsIqorGSFzRJ2efSBKq7++2KubYfdGebMzJmZM9d8P8nJue/rvs45/zkz8zvXuc593ydVhSSpL89a7AIkScNnuEtShwx3SeqQ4S5JHTLcJalDKxa7AIATTzyxxsbGFrsMSVpSdu7c+XBVrZ5s20iE+9jYGDt27FjsMiRpSUnynam2OS0jSR0y3CWpQ4a7JHXIcJekDg0U7kn2JPlqkjuT7Ghtxye5Ncm97XpVa0+SjyTZnWRXktfO5w8gSTrSTEbuv1RVp1XVhrZ+GbC9qtYD29s6wFuB9e2yGbhqWMVKkgYzl2mZjcDWtrwVOH9C+zU17nZgZZJT5vA4kqQZGjTcC/jXJDuTbG5tJ1fVgwDt+qTWvgZ4YMJt97a2p0myOcmOJDsOHDgwu+olSZMa9CCmM6tqX5KTgFuTfOMZ+maStiNOGl9VW4AtABs2bPCk8pI0RAOFe1Xta9f7k3wKOB14KMkpVfVgm3bZ37rvBdZNuPlaYN8Qa1YHxi779KTtey4/b4Erkfo07bRMkuOSPP/QMvBm4GvANmBT67YJuLEtbwPe2/aaOQN44tD0jSRpYQwycj8Z+FSSQ/3/sao+k+TLwPVJLgHuBy5s/W8GzgV2A08CFw+9aknSM5o23KvqPuDVk7Q/Apw9SXsBlw6lOknSrHiEqiR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SerQisUuQJpo7LJPT9q+5/LzFrgSaWlz5C5JHTLcJalDA4d7kqOSfCXJTW391CRfSnJvkuuSHN3aj2nru9v2sfkpXZI0lZmM3N8H3DNh/QrgyqpaDzwGXNLaLwEeq6qXAVe2fpKkBTRQuCdZC5wH/F1bD3AWcEPrshU4vy1vbOu07We3/pKkBTLo3jJ/BXwAeH5bPwF4vKoOtvW9wJq2vAZ4AKCqDiZ5ovV/eOIdJtkMbAZ48YtfPNv6NeKm2vtF0vyaduSe5G3A/qraObF5kq41wLb/b6jaUlUbqmrD6tWrBypWkjSYQUbuZwJvT3IucCzwAsZH8iuTrGij97XAvtZ/L7AO2JtkBfBC4NGhVy5JmtK0I/eq+mBVra2qMeAi4LaqejfwOeCC1m0TcGNb3tbWadtvq6ojRu6SpPkzl/3c/wB4f5LdjM+pX93arwZOaO3vBy6bW4mSpJma0ekHqurzwOfb8n3A6ZP0+SFw4RBqkyTNkkeoSlKHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SerQjE4/IC2Wqc4Lv+fy8xa4EmlpcOQuSR0y3CWpQ07LaCj8Oj1ptDhyl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjo07XeoJjkW+AJwTOt/Q1V9KMmpwLXA8cAdwHuq6qkkxwDXAK8DHgHeUVV75ql+LTC/K1VaGgYZuf8IOKuqXg2cBpyT5AzgCuDKqloPPAZc0vpfAjxWVS8Drmz9JEkLaNqRe1UV8P22+ux2KeAs4F2tfSvwJ8BVwMa2DHAD8NdJ0u5HkqY01TvDPZeft8CVLH0DzbknOSrJncB+4Fbg28DjVXWwddkLrGnLa4AHANr2J4ATJrnPzUl2JNlx4MCBuf0UkqSnGSjcq+onVXUasBY4HXjFZN3adZ5h28T73FJVG6pqw+rVqwetV5I0gBntLVNVjwOfB84AViY5NK2zFtjXlvcC6wDa9hcCjw6jWEnSYKYN9ySrk6xsy88Bfhm4B/gccEHrtgm4sS1va+u07bc53y5JC2vaD1SBU4CtSY5i/MXg+qq6KcnXgWuT/BnwFeDq1v9q4ONJdjM+Yr9oHuqWtIS5S+38G2RvmV3AayZpv4/x+ffD238IXDiU6iRJs+IRqpLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4Ocz11TmOk5qf2SX0kLxXCXnsEzvYD7Yq1RZrhLmjd+49Licc5dkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchdISWNvKl2qfRYg6k5cpekDhnuktQhp2WkWXKqQKPMkbskdciRu5YVR9t98fc5NcNdS5r/3NLknJaRpA4Z7pLUoWnDPcm6JJ9Lck+Su5O8r7Ufn+TWJPe261WtPUk+kmR3kl1JXjvfP4Qk6ekGGbkfBH6vql4BnAFcmuSVwGXA9qpaD2xv6wBvBda3y2bgqqFXLUl6RtOGe1U9WFV3tOXvAfcAa4CNwNbWbStwflveCFxT424HViY5ZeiVS5KmNKM59yRjwGuALwEnV9WDMP4CAJzUuq0BHphws72t7fD72pxkR5IdBw4cmHnlkqQpDbwrZJLnAf8E/E5VfTfJlF0naasjGqq2AFsANmzYcMR2aSH5XZ/qzUAj9yTPZjzY/6Gq/rk1P3RouqVd72/te4F1E26+Ftg3nHIlSYMYZG+ZAFcD91TVX07YtA3Y1JY3ATdOaH9v22vmDOCJQ9M3kqSFMci0zJnAe4CvJrmztf0hcDlwfZJLgPuBC9u2m4Fzgd3Ak8DFQ614CZvp0ZQefSlptqYN96r6dyafRwc4e5L+BVw6x7qWFed7JQ2bR6hKUoc8cZi65LshLXeG+xK0EHPxhqO0tDktI0kdMtwlqUOGuyR1yDl3acg8PkGjwJG7JHXIkXtHHDFqsbh31ehx5C5JHTLcJalDTstI6o5TlI7cJalLhrskdchwl6QOGe6S1CHDXZI6ZLhLUofcFXIAHn0naalx5C5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUofcz30Zc/99qV+O3CWpQ4a7JHXIaRlJy8Zy+oYmR+6S1CFH7hP0+gFjrz+XpKk5cpekDk07ck/yUeBtwP6q+pnWdjxwHTAG7AF+taoeSxLgw8C5wJPAr1fVHfNTuqSF5DvApWWQkfvfA+cc1nYZsL2q1gPb2zrAW4H17bIZuGo4ZUqSZmLacK+qLwCPHta8EdjalrcC509ov6bG3Q6sTHLKsIqVJA1mtnPuJ1fVgwDt+qTWvgZ4YEK/va3tCEk2J9mRZMeBAwdmWYYkaTLD3lsmk7TVZB2raguwBWDDhg2T9pkPzhtKWg5mO3J/6NB0S7ve39r3Ausm9FsL7Jt9eZKk2ZhtuG8DNrXlTcCNE9rfm3FnAE8cmr6RJC2cQXaF/ATwRuDEJHuBDwGXA9cnuQS4H7iwdb+Z8d0gdzO+K+TF81CzJGka04Z7Vb1zik1nT9K3gEvnWpQkaW48QlWSOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA75TUzSAllO39+pxefIXZI6ZLhLUocMd0nqkOEuSR3yA1VJy16PH3Y7cpekDhnuktQhw12SOmS4S1KH/EBV0tNM9eGilhbDXZKmsJT3onFaRpI6ZLhLUoe6nZZx3lDScubIXZI6ZLhLUocMd0nq0JKfc3duXZKO5Mhdkjq05EfukmbHd72ztxQObnLkLkkdcuQuLbKlMArU0uPIXZI65MhdkobkmT7HWOh3Yo7cJalD8zJyT3IO8GHgKODvqury+XgcSdNzr5jRsNCfrQw93JMcBfwN8CZgL/DlJNuq6uvDfiypZ37QqrmYj2mZ04HdVXVfVT0FXAtsnIfHkSRNYT6mZdYAD0xY3wu84fBOSTYDm9vq95N8cx5qma0TgYcXu4hJWNfMdFlXrhhiJU/X5fM1j4ZS1xx/ny+ZasN8hHsmaasjGqq2AFvm4fHnLMmOqtqw2HUczrpmxrpmxrpmZlTrOmQ+pmX2AusmrK8F9s3D40iSpjAf4f5lYH2SU5McDVwEbJuHx5EkTWHo0zJVdTDJbwO3ML4r5Eer6u5hP848G8npIqxrpqxrZqxrZka1LgBSdcR0uCRpifMIVUnqkOEuSR1aduGe5Jwk30yyO8llk2w/Jsl1bfuXkoxN2PaqJP+V5O4kX01y7GLXleTdSe6ccPlpktNGoK5nJ9nanqd7knxwWDXNsa6jk3ys1XVXkjcucF2/kOSOJAeTXHDYtk1J7m2XTSNU12eSPJ7kpmHWNJe6kpw24X9xV5J3jEhdL0mys/0v3p3kN4dZ14xU1bK5MP4B77eBlwJHA3cBrzysz28Bf9uWLwKua8srgF3Aq9v6CcBRi13XYX1+FrhvRJ6vdwHXtuXnAnuAsRGo61LgY235JGAn8KwFrGsMeBVwDXDBhPbjgfva9aq2vGqx62rbzgZ+BbhpWH9bQ3i+Xg6sb8svAh4EVo5AXUcDx7Tl57W/+xcN83kb9LLcRu6DnBphI7C1Ld8AnJ0kwJuBXVV1F0BVPVJVPxmBuiZ6J/CJIdU017oKOC7JCuA5wFPAd0egrlcC2wGqaj/wODCsA1Gmrauq9lTVLuCnh932LcCtVfVoVT0G3AqcMwJ1UVXbge8NqZah1FVV36qqe9vyPmA/sHoE6nqqqn7UVo9hEWdHllu4T3ZqhDVT9amqg8ATjI/SXw5Uklva27EPjEhdE72D4Yb7XOq6AfgB4yOq+4G/qKpHR6Cuu4CNSVYkORV4HU8/6G6+65qP2y7mfc/FUOpKcjrjI+Zvj0JdSdYl2dXu44r24rPgltuXdQxyaoSp+qwAfg54PfAksD3JzjaqWcy6xjcmbwCerKqvDaGeYdR1OvATxt8yrwK+mOSzVXXfItf1UeAVwA7gO8B/AgeHUNOgdc3HbRfzvudiznUlOQX4OLCpqo541zFLc6qrqh4AXpXkRcC/JLmhqh4aUm0DW24j90FOjfB/fdqUwguBR1v7v1XVw1X1JHAz8NoRqOuQixjuqH2udb0L+ExV/bhNf/wHw5v+mHVdVXWwqn63qk6rqo3ASuDeBaxrPm67mPc9F3OqK8kLgE8Df1xVt49KXYe0EfvdwM8Pqa4ZWW7hPsipEbYBh/ZUuAC4rcY/HbmF8Vfj57aw+EVgWOeon0tdJHkWcCHjc4PDNJe67gfOyrjjgDOAbyx2Xe33dxxAkjcBB2t43zUwl1Nv3AK8OcmqJKsY/4znlhGoaz7Nuq7W/1PANVX1yRGqa22S57TlVcCZwOKc8XYxPsVdzAtwLvAtxufn/qi1/Snw9rZ8LPBJYDfw38BLJ9z21xh/Jf4a8OcjVNcbgdtH6flifE+BT7bn6+vA749IXWOM/7PdA3wWeMkC1/V6xkeGPwAeAe6ecNvfaPXuBi4eobq+CBwA/qf1ecti19X+F38M3DnhctoI1PUmxvequ6tdbx7m73EmF08/IEkdWm7TMpK0LBjuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUP/C/kJUYzLN5P7AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(params[:, 1], bins=50)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([3., 1., 2., 6., 3., 2., 3., 0., 0., 3., 3., 1., 3., 1., 2., 2., 2.,\n",
" 3., 2., 1., 4., 1., 1., 1., 0., 0., 1., 1., 2., 1., 1., 2., 3., 3.,\n",
" 4., 4., 3., 5., 2., 3., 1., 2., 4., 2., 3., 0., 1., 1., 0., 1.]),\n",
" array([0.06204961, 0.06255917, 0.06306873, 0.06357828, 0.06408784,\n",
" 0.0645974 , 0.06510696, 0.06561652, 0.06612607, 0.06663563,\n",
" 0.06714519, 0.06765475, 0.06816431, 0.06867386, 0.06918342,\n",
" 0.06969298, 0.07020254, 0.07071209, 0.07122165, 0.07173121,\n",
" 0.07224077, 0.07275033, 0.07325988, 0.07376944, 0.074279 ,\n",
" 0.07478856, 0.07529812, 0.07580767, 0.07631723, 0.07682679,\n",
" 0.07733635, 0.0778459 , 0.07835546, 0.07886502, 0.07937458,\n",
" 0.07988414, 0.08039369, 0.08090325, 0.08141281, 0.08192237,\n",
" 0.08243193, 0.08294148, 0.08345104, 0.0839606 , 0.08447016,\n",
" 0.08497971, 0.08548927, 0.08599883, 0.08650839, 0.08701795,\n",
" 0.0875275 ]),\n",
" <a list of 50 Patch objects>)"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAMtklEQVR4nO3cf4xld1nH8ffTndbSpYhxrwZbrmOjNmmA0nItkBrUYpqWJeAfNbZAFaOZPwRSEiJZYvxD/9oYY8RI1ImCNQIFi42kKy0IrogpxZ26rd1uMWUzhaU1S2mU0ibW1sc/7pnd6XBn5szMPfc+e+f9Sk7m3HO/97vPM+f2kzPnRyMzkSTVdc60C5AkbcyglqTiDGpJKs6glqTiDGpJKm6ui0n37duX8/PzXUwtSTNpaWnpiczsjXqvk6Cen5/nyJEjXUwtSTMpIh5d7z1PfUhScQa1JBVnUEtScQa1JBVnUEtScQa1JBXXKqgj4qURcXtEPBwRxyPi9V0XJkkaansf9QeBuzLzhog4D7igw5okSatsGtQR8RLgDcA7ATLzWeDZbsuSJK1oc0R9CfAt4CMRcTmwBNySmU+vHhQRC8ACQL/fH3edWzZ/4NDI7csH90+4EknamTbnqOeAK4E/ycwrgKeBA2sHZeZiZg4yc9DrjXxcXZK0DW2C+iRwMjPvbV7fzjC4JUkTsGlQZ+Z/At+IiEubTW8EHuq0KknSaW3v+ngP8NHmjo8TwK92V5IkabVWQZ2ZR4FBx7VIkkbwyURJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKs6glqTiDGpJKm6uzaCIWAaeAp4HnsvMQZdFSZLOaBXUjZ/LzCc6q0SSNJKnPiSpuLZH1Al8NiIS+LPMXFw7ICIWgAWAfr8/vgolTdT8gUMjty8f3D/hSrSi7RH11Zl5JXA98K6IeMPaAZm5mJmDzBz0er2xFilJu1mroM7Mx5qfp4A7gKu6LEqSdMamQR0ReyPiwpV14Frgwa4LkyQNtTlH/cPAHRGxMv5jmXlXp1VJkk7bNKgz8wRw+QRqkSSN4O15klScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxRnUklScQS1JxbUO6ojYExH/FhF3dlmQJOmFtnJEfQtwvKtCJEmjtQrqiLgY2A/8ebflSJLWmms57g+B9wMXrjcgIhaABYB+v7/zyjTT5g8cGrl9+eD+CVcyO9b7na6n2u96o/qr1Tppmx5RR8SbgVOZubTRuMxczMxBZg56vd7YCpSk3a7NqY+rgbdExDJwG3BNRPx1p1VJkk7bNKgz8wOZeXFmzgM3Al/IzHd0XpkkCfA+akkqr+3FRAAy8zBwuJNKJEkjeUQtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScUZ1JJUnEEtScVtGtQRcX5EfCUi7o+IYxHxO5MoTJI0NNdizP8A12TmdyPiXOBLEfGZzPxyx7VJkmgR1JmZwHebl+c2S3ZZlCTpjDZH1ETEHmAJ+HHgQ5l574gxC8ACQL/f33ZB8wcOjdy+fHD/tufcyfzjrGe9udYziZq2Ylr/7m611e/LtOaf5vdit3wnW11MzMznM/PVwMXAVRHxihFjFjNzkJmDXq837joladfa0l0fmflfwGHguk6qkSR9jzZ3ffQi4qXN+ouAnwce7rowSdJQm3PULwNubc5TnwN8MjPv7LYsSdKKNnd9PABcMYFaJEkj+GSiJBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScZsGdUS8PCL+MSKOR8SxiLhlEoVJkobmWox5DnhfZt4XERcCSxHxucx8qOPaJEm0OKLOzMcz875m/SngOHBR14VJkobaHFGfFhHzwBXAvSPeWwAWAPr9/hhKa2f+wKFS81TUdW/rzb98cH+n/+40javns+l7N83/1s6m31MXWl9MjIgXA58C3puZ31n7fmYuZuYgMwe9Xm+cNUrSrtYqqCPiXIYh/dHM/NtuS5Ikrdbmro8A/gI4npl/0H1JkqTV2hxRXw3cDFwTEUeb5U0d1yVJamx6MTEzvwTEBGqRJI3gk4mSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVJxBLUnFGdSSVNymQR0RH46IUxHx4CQKkiS9UJsj6r8Eruu4DknSOjYN6sz8IvDkBGqRJI0Qmbn5oIh54M7MfMUGYxaABYB+v/+aRx99dFsFzR84tK3PTdrywf1b/szZ0tss2Or+cd/MlvX2/1b381bn2U4urIiIpcwcjHpvbBcTM3MxMweZOej1euOaVpJ2Pe/6kKTiDGpJKq7N7XkfB+4BLo2IkxHxa92XJUlaMbfZgMy8aRKFSJJG89SHJBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtScQa1JBVnUEtSca2COiKui4ivRsQjEXGg66IkSWdsGtQRsQf4EHA9cBlwU0Rc1nVhkqShNkfUVwGPZOaJzHwWuA14a7dlSZJWRGZuPCDiBuC6zPz15vXNwGsz891rxi0AC83LS4Gvjr/cidsHPDHtIibIfmfXbuoVzs5+fzQze6PemGvx4Rix7XvSPTMXgcUtFlZaRBzJzMG065gU+51du6lXmL1+25z6OAm8fNXri4HHuilHkrRWm6D+V+AnIuLHIuI84Ebg092WJUlasempj8x8LiLeDdwN7AE+nJnHOq+shpk6ldOC/c6u3dQrzFi/m15MlCRNl08mSlJxBrUkFbdrgnqzx+Aj4vsi4hPN+/dGxPyq914VEfdExLGI+PeIOL/ZfriZ82iz/NDkOtrYdvuNiLev6udoRPxfRLy6ee81Tf+PRMQfRcSoWzenoqN+Z3H/nhsRtzb78XhEfKDtnNPSUa/LzfajEXFkct1sU2bO/MLwIujXgEuA84D7gcvWjPkN4E+b9RuBTzTrc8ADwOXN6x8E9jTrh4HBtPsbZ79rxrwSOLHq9VeA1zO8t/4zwPXT7rXjfmdu/wJvA25r1i8AloH5NnPOSq/N62Vg37T7a7vsliPqNo/BvxW4tVm/HXhjc8R4LfBAZt4PkJnfzsznJ1T3du2k39VuAj4OEBEvA16Smffk8Jv+V8AvdNXAFo293+J20m8CeyNiDngR8CzwnZZzTkMXvZ51dktQXwR8Y9Xrk822kWMy8zngvxkePf8kkBFxd0TcFxHvX/O5jzR/Pv12oVMBO+l3tV/iTHBd1Myz0ZzT0kW/K2Zt/94OPA08Dnwd+P3MfLLlnNPQRa8wDPHPRsRSDP/3F6W1eYR8FrR5DH69MXPATwM/BTwDfD4iljLz88DbM/ObEXEh8CngZoZHmtO2k36Hb0a8FngmMx/cwpzT0kW/MJv79yrgeeBHgB8A/jki/qHlnNMw9l4z8wRwdWY+1lx3+FxEPJyZXxxn4eO0W46o2zwGf3pM86fS9wNPNtv/KTOfyMxngL8HrgTIzG82P58CPsbwi1HBTvpdcSMvPLo82cyz0ZzT0kW/s7p/3wbclZn/m5mngH8BBi3nnIYueiUzH2t+ngLuoM6+HWm3BHWbx+A/DfxKs34D8IXmXOzdwKsi4oLmS/AzwEMRMRcR+2B4dRl4M/AgNeykXyLiHOAXGZ4PBCAzHweeiojXNacAfhn4u27baG3s/c7w/v06cE0M7QVeBzzccs5pGHuvEbG3+SuJZvu11Nm3o037auakFuBNwH8wvIL8W8223wXe0qyfD/wN8AjDuxsuWfXZdwDHGO7M32u27QWWGN4Rcgz4IM3dIBWWHfb7s8CXR8w5aH4HXwP+mObJ1grLuPud1f0LvLjZfgx4CPjNjeassIy7V4Z3kNzfLMcq9bre4iPkklTcbjn1IUlnLYNakoozqCWpOINakoozqCWpOINakoozqCWpuP8H2ywzwL+54g4AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(params[:100, 1], bins=50)"
]
},
{
"cell_type": "code",
"execution_count": null,
"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.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
@josef-pkt
Copy link
Author

work in progress for statsmodels,
statsmodels/statsmodels#6223
statsmodels/statsmodels#6268

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment