main_test_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test over Synthetic Data - part 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will generate 1000 data points. Each data point will contain 70 predictors ($X_i = x_{i1}, x_{i2}, \\dots, x_{i70}$) and the response of the system ($y_i$). To generate the random data points $\\mathbf{X}=<X_1, X_2, \\dots, X_{1000}>$, we sample a thousand times from the binomial distribution with $n=70$ and $p=0.5$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first import the necessary modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib\n",
    "import sys\n",
    "main_spec = importlib.util.find_spec(\"main\")\n",
    "found = main_spec is not None\n",
    "project_root = \"\"\n",
    "if not found:\n",
    "    !rm -rf noisy-generalized-boolean\n",
    "    !git clone https://github.com/RohanBh/noisy-generalized-boolean\n",
    "    project_root = \"noisy-generalized-boolean/\"\n",
    "    sys.path.append(f\"{project_root}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import main\n",
    "from IPython.display import FileLink"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### NOTE\n",
    "Set the variable `load_samples` to `False` if you want to run a new MCMC (which may take a few minutes). Otherwise, the data (samples) is loaded from the previous MCMC that was run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_samples = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generate $\\mathbf{X}$, we will use a simple `numpy` trick."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 1., 0., ..., 0., 1., 0.],\n",
       "       [1., 1., 1., ..., 1., 1., 0.],\n",
       "       [1., 1., 1., ..., 1., 1., 0.],\n",
       "       ...,\n",
       "       [1., 0., 0., ..., 0., 0., 1.],\n",
       "       [0., 1., 0., ..., 1., 1., 1.],\n",
       "       [1., 1., 0., ..., 0., 0., 0.]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.random.rand(1000, 70)\n",
    "X = np.round(X)\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 70)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have $\\mathbf{X}$, it's time to generate our synthetic system's response $\\mathbf{Y}$. To do that, first we'll need to identify the noisy function of actual predictors according to which response is made. Since our statistical method works for sparse models, we say that this underlying function depends only on the following 4 predictors out of 70 predictors:\n",
    "                $$M = <x_{13}, x_{31}, x_{38}, x_{55}>$$\n",
    "\n",
    "We'll represent this function by $F: \\{0,1\\}^4 \\rightarrow \\mathbb{Z}_{3} $, where $\\mathbb{Z}_{3} = \\{0, 1, 2\\}$.\n",
    "\n",
    "We have a generated boolean function stored as a csv. Lets load it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>x55</th>\n",
       "      <th>x38</th>\n",
       "      <th>x31</th>\n",
       "      <th>x13</th>\n",
       "      <th>y0</th>\n",
       "      <th>y1</th>\n",
       "      <th>y2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7.715613e-01</td>\n",
       "      <td>4.817829e-10</td>\n",
       "      <td>2.284387e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.581223e-04</td>\n",
       "      <td>9.997419e-01</td>\n",
       "      <td>1.238610e-21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.999999e-01</td>\n",
       "      <td>7.076523e-08</td>\n",
       "      <td>2.697647e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.032021e-01</td>\n",
       "      <td>1.337596e-14</td>\n",
       "      <td>9.679786e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.925948e-01</td>\n",
       "      <td>4.670690e-03</td>\n",
       "      <td>2.734479e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.782210e-01</td>\n",
       "      <td>6.807604e-42</td>\n",
       "      <td>1.217790e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.298939e-11</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>7.649856e-39</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.986412e-01</td>\n",
       "      <td>5.024311e-07</td>\n",
       "      <td>1.358252e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.150446e-04</td>\n",
       "      <td>9.996850e-01</td>\n",
       "      <td>1.792178e-25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.755406e-26</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>1.799499e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.999991e-01</td>\n",
       "      <td>1.784450e-14</td>\n",
       "      <td>8.674651e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.474974e-01</td>\n",
       "      <td>2.493161e-01</td>\n",
       "      <td>3.031865e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.766313e-08</td>\n",
       "      <td>5.349662e-03</td>\n",
       "      <td>9.946503e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.391341e-11</td>\n",
       "      <td>9.990870e-01</td>\n",
       "      <td>9.129619e-04</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.623370e-06</td>\n",
       "      <td>3.274524e-02</td>\n",
       "      <td>9.672531e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.093945e-01</td>\n",
       "      <td>1.906055e-01</td>\n",
       "      <td>3.821038e-09</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    x55  x38  x31  x13            y0            y1            y2\n",
       "0   0.0  0.0  0.0  0.0  7.715613e-01  4.817829e-10  2.284387e-01\n",
       "1   0.0  0.0  0.0  1.0  2.581223e-04  9.997419e-01  1.238610e-21\n",
       "2   0.0  0.0  1.0  0.0  9.999999e-01  7.076523e-08  2.697647e-08\n",
       "3   0.0  0.0  1.0  1.0  9.032021e-01  1.337596e-14  9.679786e-02\n",
       "4   0.0  1.0  0.0  0.0  9.925948e-01  4.670690e-03  2.734479e-03\n",
       "5   0.0  1.0  0.0  1.0  8.782210e-01  6.807604e-42  1.217790e-01\n",
       "6   0.0  1.0  1.0  0.0  1.298939e-11  1.000000e+00  7.649856e-39\n",
       "7   0.0  1.0  1.0  1.0  9.986412e-01  5.024311e-07  1.358252e-03\n",
       "8   1.0  0.0  0.0  0.0  3.150446e-04  9.996850e-01  1.792178e-25\n",
       "9   1.0  0.0  0.0  1.0  2.755406e-26  1.000000e+00  1.799499e-17\n",
       "10  1.0  0.0  1.0  0.0  9.999991e-01  1.784450e-14  8.674651e-07\n",
       "11  1.0  0.0  1.0  1.0  4.474974e-01  2.493161e-01  3.031865e-01\n",
       "12  1.0  1.0  0.0  0.0  1.766313e-08  5.349662e-03  9.946503e-01\n",
       "13  1.0  1.0  0.0  1.0  4.391341e-11  9.990870e-01  9.129619e-04\n",
       "14  1.0  1.0  1.0  0.0  1.623370e-06  3.274524e-02  9.672531e-01\n",
       "15  1.0  1.0  1.0  1.0  8.093945e-01  1.906055e-01  3.821038e-09"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(f\"{project_root}data/bool_func.csv\", index_col=0)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each row of this boolean function is sampled from a dirichlet distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define a utility function that generates the response for a row in $\\mathbf{X}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bool_func(xi):\n",
    "    \"\"\"Takes an index/row of truth table and returns the response \n",
    "    according to the boolean function\"\"\"\n",
    "    if isinstance(xi, int):\n",
    "        weight_arr = df.iloc[xi, 4:].values\n",
    "    else:\n",
    "        bool_arr = np.all(df.iloc[:, :4].values == xi, axis=1)\n",
    "        weight_arr = df.iloc[bool_arr, 4:].iloc[0].values\n",
    "    return np.random.choice(3, p=weight_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `bool_func` we can start generating responses of our synthetic system. The response depends only on **actual predictors** of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 4)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_actual = X[:, [55, 38, 31, 13]]\n",
    "X_actual.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we generate the response of the synthetic system $\\mathbf{Y}$, for all $1000$ rows in $\\mathbf{X}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000,)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = np.apply_along_axis(bool_func, axis=1, arr=X_actual)\n",
    "Y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have all we need. Let's see if our statistical method can predict which predictors does the response depends upon.\n",
    "\n",
    "We call the `metropolis` function from the module `main`. This function uses a **Metropolis-Hastings** MCMC algorithm. This algorithm runs a markov chain over the model space $\\mathcal{M}$ and generates samples of possible models with the posterior probability $P(M \\mid \\mathbf{Y}, \\mathbf{X})$. This behavior is true only after the markov chain has reached a steady state.\n",
    "\n",
    "The `metropolis` function returns the generated samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "********** Running MCMC **********\n",
      "500: 500 iterations took 42.59 seconds                                         \n",
      "1000: 500 iterations took 51.15 seconds                                         \n",
      "1500: 500 iterations took 25.43 seconds                                         \n",
      "2000: 500 iterations took 25.37 seconds                                         \n",
      "2500: 500 iterations took 25.12 seconds                                         \n",
      "3000: 500 iterations took 23.04 seconds                                         \n",
      "3500: 500 iterations took 15.94 seconds                                         \n",
      "4000: 500 iterations took 15.97 seconds                                         \n",
      "4500: 500 iterations took 15.50 seconds                                         \n",
      "5000: 500 iterations took 15.15 seconds                                         \n",
      "5500: 500 iterations took 15.06 seconds                                         \n",
      "6000: 500 iterations took 15.09 seconds                                         \n",
      "6500: 500 iterations took 15.18 seconds                                         \n",
      "7000: 500 iterations took 15.27 seconds                                         \n",
      "7500: 500 iterations took 15.23 seconds                                         \n",
      "********** loop took 346.26 seconds **********                               \n"
     ]
    }
   ],
   "source": [
    "samples = []\n",
    "if not load_samples:\n",
    "    samples = main.metropolis(X, Y, q=3, iterations=8000, show_progress=True)\n",
    "    pd.DataFrame(samples).to_csv(f\"{project_root}data/mcmc-1.csv\")\n",
    "else:\n",
    "    samples = pd.read_csv(f\"{project_root}data/mcmc-1.csv\", index_col=0).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a href='data/mcmc-1.csv' target='_blank'>data/mcmc-1.csv</a><br>"
      ],
      "text/plain": [
       "/home/rohan/Projects/PythonProjects/generalized-noisy-boolean/data/mcmc-1.csv"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FileLink(f'{project_root}data/mcmc-1.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Click on the link above to download the generated samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ioff() # turn off interactive mode for matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8001, 70)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = np.array(samples)\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHiNJREFUeJzt3XucHGWd7/HPt2cmFyDkQgaIhjCgKKJICEO4iQdREdGjyKILrCysnldcr3jZdUGPiq+jHnQRleMF2EXBPYgiF+GAclkEWUTBBEISSDCgQQIhGYHAEMhlJr/zRz2dTE96Mt2T6a4K/X2/XvPq6urqql/SnXznqeeppxQRmJmZlZXyLsDMzIrFwWBmZhUcDGZmVsHBYGZmFRwMZmZWwcFgZmYVHAxmZlbBwWBmZhUcDGZmVqG9WQeSNA64AxibjntlRHxJ0iXAfwOeTZueHhHzt7avqVOnRldXVwOrNTN76Zk3b95fI6JzuO2aFgzAOuDoiHheUgdwp6Rfpdf+OSKurHVHXV1dzJ07tyFFmpm9VEl6tJbtmhYMkU3K9Hx62pF+PFGTmVnBNLWPQVKbpPnAKuCWiLg7vfRVSQskfUvS2GbWZGZmlZoaDBHRHxEzgenAbEmvA84C9gUOBqYA/1LtvZLmSJoraW5PT0/TajYzazW5jEqKiNXA7cCxEbEiMuuAHwGzh3jPRRHRHRHdnZ3D9p2YmdkINS0YJHVKmpSWxwNvAZZImpbWCTgeWNSsmszMbEvNHJU0DbhUUhtZIF0REddL+rWkTkDAfOAfm1iTmZkN0sxRSQuAA6usP7pZNZiZ2fBa6srntRv62efzv+Rz1yzEtzQ1M6uupYLh6HNvZ0N/8JO7/8Ltf/TIJjOzaloqGJ54du2m5efX9uVYiZlZcbVUMJiZ2fAcDGZmVqFlg8Fdz2Zm1bVsMJiZWXUtGwy9azfkXYKZWSG1bDC0l5R3CWZmhdSyweDr28zMqmupYGgb0ErY6GAwM6uqpYJh4Omj8LgkM7OqWioY3GIwMxteSwVDSQM6nN3JYGZWVUsFw4kHTd+07BaDmVl1LRUMX3jnfpuWPe22mVl1LRUM7mMwMxteSwXDQM4FM7PqWjcYfCrJzKyqFg6GvCswMyumlg2Gq+5dnncJZmaF1LLBsOTJ3rxLMDMrpJYLhoVnH5N3CWZmhdZywVBx9bOZmW2hacEgaZykeyTdL+kBSV9O6/eSdLekpZJ+JmlMY+to5N7NzLZ/zWwxrAOOjogDgJnAsZIOBb4OfCsi9gGeAT7YyCLcYjAz27qmBUNknk9PO9JPAEcDV6b1lwLHN6smMzPbUlP7GCS1SZoPrAJuAR4BVkdEX9pkOfDyRtbgFoOZ2dY1NRgioj8iZgLTgdnAa6ptVu29kuZImitpbk9Pz4hrcC6YmW1dLqOSImI1cDtwKDBJUnt6aTrwxBDvuSgiuiOiu7Ozc8THdi6YmW1dM0cldUqalJbHA28BFgO3ASemzU4Drm1kHT6VZGa2de3DbzJqpgGXSmojC6QrIuJ6SQ8CP5X0FeA+4OJGFuFcMDPbuqYFQ0QsAA6ssv5PZP0NTSEng5nZVrXclc9mZrZ1DgYzM6vQksHwjv2n8cpdd8q7DDOzQmrJYJBgo+/UY2ZWVUsGQ0nyHdzMzIbQosHgFoOZ2VBaNBjkYDAzG0IzL3ArjF8tepIXN/SzZl0fO45tyb8CM7MhtWSL4cUN/QA8vvrFnCsxMyuelgwGMzMbmoPBzMwqOBjMzKyCg8HMzCo4GMzMrEJLB4MvZTAz21JLB4OZmW2ppYPB9+wxM9tSSweDTyWZmW2ppYPB8yWZmW3JwWBmZhVaOhicC2ZmW2rJYDjjzfsAbjGYmVXTksEwc49JAPRvdDCYmQ3WksFQKmXjVJ0LZmZbalowSNpD0m2SFkt6QNIZaf3Zkh6XND/9HNfoWlIuED6VZGa2hWbevqwP+ExE3CtpAjBP0i3ptW9FxLnNKqQktxjMzIbStGCIiBXAirTcK2kx8PJmHX+g8hXP6/s25nF4M7NCy6WPQVIXcCBwd1r1MUkLJP1Q0uRGH7/cYvjUFfMbfSgzs+1O3cEgaUdJbSM9oKSdgKuAT0bEc8APgFcAM8laFN8c4n1zJM2VNLenp2ekhwc2X7/Q07tum/ZjZvZSNGwwSCpJOkXSDZJWAUuAFakD+V8l7VPrwSR1kIXCZRFxNUBErIyI/ojYCPwbMLvaeyPioojojojuzs7OWg9ZlTudzcyGVkuL4Tay3+jPAnaPiD0iYlfgSOD3wDmS3j/cTiQJuBhYHBHnDVg/bcBm7wEW1VH/iPQ7GMzMhlRL5/NbImLD4JUR8TTZb/9XpZbAcI4ATgUWSiqf3P8ccLKkmUAAy4AP1VL4tvBoJDOzoQ0bDNVCYYTb3AlUuwPCL4d772jzVBhmZkMbNhgkzahxX6tTZ3LhuY/BzGxotZxKurSGbQK4BPjxNlXTJAfukY2IHd8x4sFVZmYvWbUEw/8CFkbEto0RLZDJO45hdtcU2kq+t6eZ2WC1BMN/AqskbSQbMbQAWJgeH4yI7fJiAMl9DWZm1dQyXPUTwBPA+cBXgIeAg4DzgEcbV1pjlSTfqMfMrIphgyEivks21DSAbwMbgDMi4k0RsXuD62uYUsktBjOzamqaEiMiXoyIrwNHAa8E7pF0SCMLa7SS5GAwM6uiluGqRwKvAfZNj7sCvcAujS2tsST5Qjczsypq6Xz+DXA/cDlwfkQsa2hFTVKSr2cwM6umlmD4MLA/8A7gM5KeIhuVtBBYFBG/aGB9DVNyi8HMrKpapsS4cOBzSdOB15OFxd8A22kwuPPZzKyauu/gFhHLgeXkMMfRaHIfg5lZdS05VxK4j8HMbCi1zpUUVJ8ZtWy7misJ4Ok161nyZG/eZZiZFU4tfQxvakYhzfaHZc8AsHZDP+M8mZ6Z2Sa13Npz2JnmatmmqPrc0WBmVqGmW3tK+vjgvgZJYyQdLelS4LTGlNd4HplkZlaplj6GY4EPAJdL2gtYDYwD2oCbgW9FxPytvL/QYmPeFZiZFUstfQxrge8D30/3dp4KvBgRqxtdXDP0u8VgZlahrusY0r2dVzSollz4VJKZWaWaZld9KdvozmczswotGwzvPWg6gK9+NjMbpOZgUOb9kr6Yns+QNLtxpTVWd9dkwH0MZmaD1dNi+D5wGHByet4LfG/UK2qSUrr0wqeSzMwq1RMMh0TER4G1ABHxDDCm1jdL2kPSbZIWS3pA0hlp/RRJt0hamh4n1/UnGKFyMLjBYGZWqZ5g2CCpjWxeJCR1AvVcBdAHfCYiXgMcCnxU0n7AmcCtEbEPcGt63nCl9Cf3qSQzs0r1BMP5wDXArpK+CtwJ/O9a3xwRKyLi3rTcCywGXg68m2yiPtLj8XXUNGKbTiU5GMzMKtR8HUNEXCZpHvBmsplWj4+IxSM5qKQu4EDgbmC3iFiRjrFC0q5DvGcOMAdgxoxaZwIfmvsYzMyqq2dU0q3A3hHxvYj4bkQslnRRvQeUtBNwFfDJeu7fEBEXRUR3RHR3dnbWe9gttJXKLYZt3pWZ2UtKPaeS9gL+RdKXBqzrrudgaUqNq4DLIuLqtHqlpGnp9WnAqnr2OVIpF+h3MpiZVagnGFaTnUbaTdL/kzSxngOlqbkvBhZHxHkDXrqOzbOzngZcW89+R8p9DGZm1dUzV5Iiog/4iKTTyTqf6xlaegRwKrBQUnk21s8B5wBXSPog8BfgvXXsc8Q8XNXMrLp6guGC8kJEXCJpIfDRWt8cEXcy9O1B31xHHaPCw1XNzKqr+VRSRFw46Pm8iPjA6JfUHOUWwxevXZRzJWZmxVLLrT3vTI+9kp5LP73l540vsTF2nzgOgAXLn825EjOzYqnlRj1vSI8TGl9O80ybOH7TckSwHd+22sxsVNVzHcN7JU1Iy/9T0tWSDmxcaY3V0bY5CDxk1cxss3qGq34hInolvQF4G9n0FRcM857CKl/gBu6ANjMbqJ5g6E+P7wB+EBHXUsfsqkXTUdr8R99Yz1SAZmYvcfUEw+OSLgTeB/xS0tg6318oJbcYzMyqquc/9vcBNwHHRsRqYArwzw2pqsn+sOzpvEswMyuMeq5jeCEiro6Ipen5ioi4uXGlNc9tS5oyPZOZ2XZhuz0VNJo8KsnMbDMHA55Iz8xsIAcD0NfvYDAzKxv2ymdJvaT7PLN5ErxIyxEROzeotqbxqCQzs81qmRLjJTUVRjXzHn0m7xLMzAqjnikxJOn9kr6Qnu8haXbjSmueR596Ie8SzMwKo54+hu8DhwGnpOfPA98b9Yqa6KJTD8q7BDOzwqnnRj2HRMQsSfcBRMQzkrbbKTEA2ts8o6qZ2WD1tBg2SGojdURL6gS261mGSp5q28xsC/UEw/nANcCukr5Kds/nrzWkqiYZOMOqmZllaj6VFBGXSZpHdn9mAcdHxOKGVdYEbjGYmW2pnj4GImIJsKRBtTTdmHZf32dmNlgtF7h9emuvR8R5o1dOc3XvORmAKTtu133oZmajqpYWQ/kCt1cDBwPXpef/HbijEUU1iyQOmD6RSTs4GMzMymq58vnLAJJuBmZFRG96fjbw84ZW1wSlkjyJnpnZAPWcZJ8BrB/wfD3QVeubJf1Q0ipJiwasO1vS45Lmp5/j6qhnVJTkYDAzG6iezuf/AO6RdE16fjxwaR3vvwT4LvDjQeu/FRHn1rGfUdUm+Z7PZmYD1DNc9auSfgUcSXaR2z9ExH11vP8OSV11V9hga/v6WbD8WR5e1cufetYwdcJYZs2YzINPPMe+u0+gVBLzHn2a59b2MbatxOumT2TncR0ArO/byMLHn2XWjEk8t7aP5c+8wGtfNjHnP5GZ2bapa7gq0E92tXMwelc9f0zS3wNzgc9ERNWpTiXNAeYAzJgxY5QODQuWPwvAW867Ix0HfvyB2Zx68T185q2v4uRDZvA3P/jdpu3fPfNlfOekAwG4+M4/8/Ubl3D1Rw7nwt88wk0PrOSRrx3nC+fMbLtWz+yqZwCXAVOBXYH/K+nj23j8HwCvAGYCK4BvDrVhRFwUEd0R0d3Z2bmNhx1aBDyy6nkA7l/+LGvW9VW8PnfZ5ty69y/Zck/vOm56YCUAfT4vZWbbuXpaDB8km0hvDYCkrwO/A/7PSA8eESvLy5L+Dbh+pPsaTeUbukUEfXXeD7qvPxhbbzvMzKxA6hmVJLJTSWX9bL6j24hImjbg6XuARUNt20wbUxgE0F9vMNS5vZlZ0dTzu+2PgLvTqCQB7wZ+WOubJV0OHAVMlbQc+BJwlKSZZP8HLwM+VEc9DVO+1WdEsKG/vlNDfXVub2ZWNPWMSjpP0u3AG9Kq0yNifh3vP7nK6otrfX8zlVsJG6P+FkO925uZFY1imIu7JF03eFV6DICIeFcD6tqq7u7umDt37qjsq+vMG7ZYt9vOY1n53DoADthjEvc/trri9SP3mQrAfy396xbv3XXCWF65605MHN/BOSe8nok7ZENbb3toFT+75zG+dsL+npvJzHIhaV5EdA+3XS0thsOAx4DLgbvZxn6Fojn3vQfwTz+/H4DTD+9i6ape1qzrZ+Vz65g+eTwlwRGv3IX2Uonf/LGHV3TuyPNppNL4jjZe3NBfsb9VvetY1ZuFyqmH7cnhr8hC5Bs3PsTiFc9xyiEzeOOrGjeqysxsW9USDLsDbwVOJrvf8w3A5RHxQCMLa5YTD5rOiQdNH9F7H3v6BY78xm1Dvj7wtNK6FCA+1WRmRTfsqKSI6I+IGyPiNOBQ4GHg9lG4hmG7N9yFbH39W4ZAvZ3ZZmbNVlPns6SxwDvIWg1dZLf5vLpxZW0f2tuGCYYqrQO3GMys6Gq5Uc+lwOuAXwFfjohCXGtQBO2lrTe4+qtcBb3BwWBmBVdLi+FUYA3wKuAT2nyfZAERETs3qLbCG+5U0iM9a3j0qTUA9DyfdUivfHbtpnUAk3YYw8TxHY0r0sysTsMOVy2i0Ryuui3Wbuhn3y/cuE372HlcO/d98RhPvGdmDTeaw1VtCOM62vjpnEMZ19HGM2vW07uuj7se/ms20umCbEbW8953AFfOW85djzzFzuPaOftdr930/luXrOKGBStY37eR8WPa8vpjmJlVcDBso0P33qXi+bsOeFnF8xNmTef+x1Zz1yNPcfrhXZwwa/PQ2KeeX88NC1akGVkdDGZWDPVMomcjVEqnidoGdVaXRzV5pJKZFYmDoYkGD29tT4Gxocr1DmZmeXEwNFH7oA7m9rbsr98tBjMrEgdDEw0eeVR+/vjqF7n0rmUOCDMrBHc+N8g7Xz+NjtQiOLhrClfOXc6rd59Qsc3kHbJZVj99xXwefeoFOtpKnHLI6N3P2sxsJBwMDfLdU2ZtWj5u/2kct/+0LbZ57cuyawOfen49AC+s79tiGzOzZvOppByV+xzW92VTZwy4qtzMLDcOhhyVO583pDmVfPGzmRWBgyFH5c7n8qwkJbcYzKwAHAw56hh0XUPJTQYzKwAHQ44GT9vd5haDmRWAgyFHY9pLfPU9r2P2XlMA2H3i2JwrMjNzMOTu7w7Zky+8Yz8AqtzXx8ys6RwMBVA+g7RxO7w3hpm99DQtGCT9UNIqSYsGrJsi6RZJS9Pj5GbVUyTl0UieEcPMiqCZLYZLgGMHrTsTuDUi9gFuTc9bTrkPenu8m56ZvfQ0LRgi4g7g6UGr3w1cmpYvBY5vVj1FUm4x9DsYzKwA8u5j2C0iVgCkx12H2lDSHElzJc3t6elpWoHN4FNJZlYkeQdDzSLioojojojuzs7OvMsZVeXr2nwqycyKIO9gWClpGkB6XJVzPbnY3GJwMJhZ/vIOhuuA09LyacC1OdaSm03B4OsYzKwAmjlc9XLgd8CrJS2X9EHgHOCtkpYCb03PW075OgZ3PptZETTtRj0RcfIQL725WTUU1eZZVh0MZpa/vE8lGR6VZGbF4mAogJKnxDCzAnEwFIDcYjCzAnEwFMCmFoOTwcwKwMFQAOXOZ59KMrMicDAUgE8lmVmROBgKoHwq6bGnX6B37YZ8izGzludgKICOtuxjuOSuZbz3gt/lXI2ZtToHQwGM62hjr6k7ArDkyd6cqzGzVudgKIjXTJuQdwlmZoCDoTDaSv4ozKwY/L9RQXSUe6DNzHLmYCiINgeDmRWEg6Eg2ts2B8NZVy/gynnLc6zGzFqZg6EgDtpzyqbla+57nG/e/FCO1ZhZK2va/Rhs6048aDonHjQdgLOuXsgtD67MuSIza1VuMRRQe0n0+z6fZpYTB0MBtbeJvn5PnGRm+XAwFFB7SfR5Rj0zy4mDoYDa20ps6PepJDPLh4OhgMa0lejbGHznP5fmXYqZtSAHQwGdNHsPAJau8oR6ZtZ8DoYCmjZxPK/abSd3QJtZLhwMBdVWKrkD2sxyUYgL3CQtA3qBfqAvIrrzrSh/HW2+lsHM8lGIYEjeFBF/zbuIomjzkFUzy4lPJRVUR6nEfy39K30etmpmTVaUYAjgZknzJM2ptoGkOZLmSprb09PT5PKab2NkrYX7HludcyVm1mqKEgxHRMQs4O3ARyW9cfAGEXFRRHRHRHdnZ2fzK2yyzx67LwBrN/TnXImZtZpCBENEPJEeVwHXALPzrSh/Y9qzj8b9DGbWbLkHg6QdJU0oLwPHAIvyrSp/7emObr6WwcyarQijknYDrpEEWT0/iYgb8y0pf+U7unnIqpk1W+7BEBF/Ag7Iu46iKbcYNrjFYGZNlnswWHVtpews31dueJDzb/VkemaW+doJ+3Nw15ThN9wGDoaC2mPyeN5/6AyeXrM+71LMrEDGd7Q1/BgOhoJqbyvxleP3z7sMM2tBuY9KMjOzYnEwmJlZBQeDmZlVcDCYmVkFB4OZmVVwMJiZWQUHg5mZVXAwmJlZBUVsf3PxSOoBHh3h26cCRbyFqOuqj+uqj+uqT1Hrgm2rbc+IGPaGNttlMGwLSXMjojvvOgZzXfVxXfVxXfUpal3QnNp8KsnMzCo4GMzMrEIrBsNFeRcwBNdVH9dVH9dVn6LWBU2oreX6GMzMbOtascVgZmZb0VLBIOlYSQ9JeljSmU043g8lrZK0aMC6KZJukbQ0PU5O6yXp/FTbAkmzBrzntLT9UkmnbWNNe0i6TdJiSQ9IOqMIdaX9jZN0j6T7U21fTuv3knR3Os7PJI1J68em5w+n17sG7OustP4hSW8bhdraJN0n6fqi1JT2uUzSQknzJc1N64rwWU6SdKWkJem7dljedUl6dfp7Kv88J+mTedeV9vep9J1fJOny9G8hv+9YRLTED9AGPALsDYwB7gf2a/Ax3wjMAhYNWPcN4My0fCbw9bR8HPArQMChwN1p/RTgT+lxclqevA01TQNmpeUJwB+B/fKuK+1TwE5puQO4Ox3zCuCktP4C4MNp+SPABWn5JOBnaXm/9PmOBfZKn3vbNtb2aeAnwPXpee41pf0uA6YOWleEz/JS4H+k5THApCLUNaC+NuBJYM+86wJeDvwZGD/gu3V6nt+xbf4L3l5+gMOAmwY8Pws4qwnH7aIyGB4CpqXlacBDaflC4OTB2wEnAxcOWF+x3SjUdy3w1gLWtQNwL3AI2cU87YM/R+Am4LC03J620+DPduB2I6xlOnArcDRwfTpGrjUN2M8ytgyGXD9LYGey/+hUpLoG1XIM8Nsi1EUWDI+RBU17+o69Lc/vWCudSir/5ZctT+uabbeIWAGQHndN64eqr2F1pybogWS/mReirnTKZj6wCriF7Lee1RHRV+U4m2pIrz8L7NKA2r4NfBbYmJ7vUoCaygK4WdI8SXPSurw/y72BHuBH6fTbv0vasQB1DXQScHlazrWuiHgcOBf4C7CC7Dszjxy/Y60UDKqyrkhDsoaqryF1S9oJuAr4ZEQ8V5S6IqI/ImaS/ZY+G3jNVo7T8NokvRNYFRHzBq7Os6ZBjoiIWcDbgY9KeuNWtm1Wbe1kp1B/EBEHAmvITtHkXVd2sOxc/buAnw+3aTPqSn0a7yY7/fMyYEeyz3OoYzS8rlYKhuXAHgOeTweeyKGOlZKmAaTHVWn9UPWNet2SOshC4bKIuLoodQ0UEauB28nO7U6S1F7lOJtqSK9PBJ4e5dqOAN4laRnwU7LTSd/OuaZNIuKJ9LgKuIYsTPP+LJcDyyPi7vT8SrKgyLuusrcD90bEyvQ877reAvw5InoiYgNwNXA4OX7HWikY/gDsk3r6x5A1Ja/LoY7rgPIohtPIzvGX1/99GglxKPBsatbeBBwjaXL6zeKYtG5EJAm4GFgcEecVpa5UW6ekSWl5PNk/mMXAbcCJQ9RWrvlE4NeRnVy9Djgpjd7YC9gHuGckNUXEWRExPSK6yL4zv46Iv8uzpjJJO0qaUF4m+wwWkfNnGRFPAo9JenVa9WbgwbzrGuBkNp9GKh8/z7r+AhwqaYf077P895Xfd2w0OnK2lx+yUQZ/JDtv/fkmHO9ysnOGG8jS/INk5wJvBZamxylpWwHfS7UtBLoH7OcDwMPp5x+2saY3kDUvFwDz089xedeV9vd64L5U2yLgi2n93ukL/jBZ839sWj8uPX84vb73gH19PtX8EPD2Ufo8j2LzqKTca0o13J9+Hih/pwvyWc4E5qbP8hdko3eKUNcOwFPAxAHrilDXl4El6Xv/H2Qji3L7jvnKZzMzq9BKp5LMzKwGDgYzM6vgYDAzswoOBjMzq+BgMDOzCg4Ga0mSnk+PXZJOGeV9f27Q87tGc/9mjeZgsFbXBdQVDJLahtmkIhgi4vA6azLLlYPBWt05wJHK5uf/VJrE718l/SHNwf8hAElHKbuPxU/ILnZC0i/S5HUPlCewk3QOMD7t77K0rtw6Udr3ImX3UPjbAfu+XZvvX3BZugIWSedIejDVcm7T/3asJbUPv4nZS9qZwD9FxDsB0n/wz0bEwZLGAr+VdHPadjbwuoj4c3r+gYh4Ok3f8QdJV0XEmZI+FtlEgIOdQHZF8AHA1PSeO9JrBwKvJZvb5rfAEZIeBN4D7BsRUZ4uxKzR3GIwq3QM2fw488mmI9+FbM4ZgHsGhALAJyTdD/yebPKyfdi6NwCXRzaD7ErgN8DBA/a9PCI2kk1T0gU8B6wF/l3SCcAL2/ynM6uBg8GskoCPR8TM9LNXRJRbDGs2bSQdRTbJ32ERcQDZHE/jatj3UNYNWO4nu0FLH1kr5SrgeODGuv4kZiPkYLBW10t2i9Oym4APp6nJkfSqNHPpYBOBZyLiBUn7kk0PXrah/P5B7gD+NvVjdJLd+nXI2S+V3TNjYkT8Evgk2Wkos4ZzH4O1ugVAXzoldAnwHbLTOPemDuAest/WB7sR+EdJC8hmsvz9gNcuAhZIujeyKbrLriG7ReP9ZDPcfjYinkzBUs0E4FpJ48haG58a2R/RrD6eXdXMzCr4VJKZmVVwMJiZWQUHg5mZVXAwmJlZBQeDmZlVcDCYmVkFB4OZmVVwMJiZWYX/D/uouoi/hCacAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(data.shape[0]), np.count_nonzero(data, axis=1))\n",
    "plt.xlabel(\"Iterations\")\n",
    "plt.ylabel(r\"Model size ($\\vert M \\vert$)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the markov chain enters steady state after around 2500 iterations. We discard the samples before that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5501, 70)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = data[2500:]\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next plot gives the posterior probability that predictor $i$ is a part of the model given $Y$ and $X$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "     \\begin{align*}\n",
    "     P(j \\in M \\mid Y, X) &= \\sum_{M \\in \\mathcal{M}} P(j \\in M \\mid \\mathbf{Y}, \\mathbf{X}, M) P(M \\mid \\mathbf{Y}, \\mathbf{X}) \\\\\n",
    "         &= \\hat{\\mathbb{E}}[I_M(j)]_{P(M|\\mathbf{Y},\\mathbf{X})} \\\\\n",
    "         &= \\frac{1}{T} \\sum_{t=1}^T I_{M^{(t)}} (j)\n",
    "     \\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEOCAYAAAB8aOvdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFR9JREFUeJzt3XvUXXV95/H3Ry4KiFJMGBkSDGqKUOWiGdRSW0TsCrQD49R2oGjVWrOmSqsL2g6MlTrMmo7VDjousZWFDh2nAyJVzKKxqBR1lhYkEC5JMDVDQZ6BMVHwjnLxO3/sHefw5PklJ5dzSfJ+rXXW2ft3fs/e35yc5PPs3z77t1NVSJI0lydNugBJ0vQyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlq2nvSBeyoefPm1aJFiyZdhiTtUm655ZZvVtX8rfXb5UNi0aJFrFy5ctJlSNIuJcm9w/RzuEmS1GRISJKaDAlJUpMhIUlqMiQkSU1jC4kkH0myIcnqxutJ8v4k65PckeSF46pNkjS3cR5JXA4s3cLrpwKL+8cy4C/GUJMkaQvGFhJV9UXgwS10OQP479W5ETgoyaHjqU6SNJdpOidxGHDfwPpM3yZJmpBpuuI6c7TVnB2TZXRDUhx++OGjrGmXt+j8v/3p8j3v+pUJVrJrGXzfYLzvnX9nw5nk39GeZJqOJGaAhQPrC4D75+pYVZdW1ZKqWjJ//lanHpEkbadpConlwG/133J6CfCdqnpg0kVJ0p5sbMNNSa4ATgLmJZkB/gTYB6Cq/hJYAZwGrAd+CLxhXLVJkuY2tpCoqrO28noBbxlTOZKkIUzTcJMkacoYEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLUZEhIkpoMCUlSkyEhSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQkSU2GhCSpyZCQJDUZEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqWmsIZFkaZJ1SdYnOX+O1w9PckOSVUnuSHLaOOuTJD3R2EIiyV7AJcCpwNHAWUmOntXtj4Grqup44Ezgg+OqT5K0uXEeSZwArK+qu6vqEeBK4IxZfQp4Wr/8dOD+MdYnSZpl7zHu6zDgvoH1GeDFs/q8E/hMkt8DDgBOGU9pkqS5jPNIInO01az1s4DLq2oBcBrw0SSb1ZhkWZKVSVZu3LhxBKVKkmC8ITEDLBxYX8Dmw0lvBK4CqKp/AJ4CzJu9oaq6tKqWVNWS+fPnj6hcSdI4Q+JmYHGSI5LsS3dievmsPl8HXgGQ5Ci6kPBQQZImZGwhUVWPAecA1wF30X2LaU2Si5Kc3nc7D3hTktuBK4DXV9XsISlJ0piM88Q1VbUCWDGr7cKB5bXAieOsSZLU5hXXkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQkSU2GhCSpyZCQJDUZEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLUZEhIkpoMCUlS0zaHRJIDkuw1imIkSdNlqyGR5ElJfjPJ3ybZAHwVeCDJmiTvSbJ49GVKkiZhmCOJG4DnABcAz6yqhVV1CPAy4EbgXUleM8IaJUkTsvcQfd5UVf84u7GqHkxyb1X9WpJ9RlCbJGnChjmSuDXJvxs8D5HkwCQfAP4BoKoeHVWBkqTJGSYkZoD/DNyW5GVJfp3uvMSbgS+PsjhJ0mQNExIvAN5Jd17i88CVwL7A71TVL42sMknSxG01JPqhpE8B64H0j1XAZ0dbmiRp0ob5CuyfAzcDRwJ/Dvw1cAqwNsl5oy1PkjRJwww3nQt8BXhRVf1RVb0WeCXwAPDuURYnSZqsYb4Cu6yqLhtsqKrrkzwf+PejKUuSNA22GhKzA2Kg/RG6E9qSpN2UE/xJkpq2OySSHJrkydv4M0uTrEuyPsn5jT6/kWRtPzfU/9ze+iRJO26YcxItHwWek+RvquoPtta5v2L7ErqT3jPAzUmWV9XagT6L6eaIOrGqHkpyyA7UJ0naQdsdElV1SpIARw/5IycA66vqboAkVwJnAGsH+rwJuKSqHur3sWF765Mk7bhhrpPYr/VaddYMua/DgPsG1mf6tkE/C/xski8luTHJ0iG3LUkagWHOSXw1yb/eCfvKHG01a31vYDFwEnAWcFmSgzbbULIsycokKzdu3LgTSpMkzWWYkHgm8PEkn07y3B3Y1wywcGB9AXD/HH0+VVWPVtU/AevoQuMJqurSqlpSVUvmz5+/AyVJkrZkmHMSPwdcDPwqcGeSv2Rg2KiqLh5yXzcDi5McAfwf4EzgN2f1uYbuCOLyJPPohp/uHnL7kqSdbJiL6dYDpye5BPhd4K2DL9MFyFZV1WNJzgGuA/YCPlJVa5JcBKysquX9a7+cZC3wOPCHVfWtbfoTSZJ2mq2GRJLjgPfR3a70e3ST/H19e3ZWVSuAFbPaLhxYLrq5os7dnu1LknauYYabVtKdu/gYcG5VPTDakiRJ02KYkFgHnFNVN4y6GEnSdBkmJI6tqsdGXokkaeoMc2c6A0KS9lDOAitJajIkJElNw3wF9grgWVvqQvft1Z/faVVJkqbCMBfTnTWOQiRJ08fhJklSkyEhSWoa5n4SJydxqlVJ2gMNczHd54ANSX4CrAbuAO7sn9dW1Y9HWJ8kaYKGCYnfB34buAr4MnAk8CLg9cBRdPebkCTthoa54voDwIl004K/D3gUeGtVvbyqDAhJ2o0NdeK6qh6uqj+ju63oc4GvJHnxKAuTJE3eMBfTvYxuWOl5/fMhdPeVeMZoS5MkTdow5yS+ANwOXAG8v6ruGWlFkqSpMUxI/C7wAuBXgPOSfIvu2013Aqur6poR1idJmqBhpuX40OB6kgXAMXTB8WuAISFJu6lhjiSeoKpmgBlm3atakrT7cVoOSVKTISFJahoqJJJcM7B87VztkqTdz7BHEocNsSxJ2s0MGxI1xLIkaTfjOQlJUtPQIZFk/02L/fpTR1KRJGlqDBsSlwIfSBKg+uf39+2SpN3UsLPAXkZ3w6GbgWcDtwB39u2SpN3U0FdcV9XFSd4HzAc2VtVPRleWJGkabNO0HH0wfGNEtUiSpozfbpIkNW01JJL8af/8i0kO2pGdJVmaZF2S9UnO30K/VyepJEt2ZH+SpB0zzHDTVf3zm4Bj+q/C3kV3I6LbgZuq6r6tbSTJXsAlwCvpZpG9Ocnyqlo7q9+BwO8DNw39p5AkjcQww023A1TVa6vqWOBo4B3A14CXAp9M8q4htnMCsL6q7q6qR4ArgTPm6PcfgXcDPxpim5KkERomJG5I8ntJDgeoqker6na6/+RX0B1V/M4Q2zkMGDzimGHW3E9JjgcWVtW1SJImbpiQWAo8DlyR5IEka5PcTXckcSbwX4DDh9hO5mj76dxPSZ4EvBc4b6sbSpYlWZlk5caNG4fYtSRpewxz+9IfAR8EPphkH2Ae8HBVfXsb9zUDLBxYXwDcP7B+IPB84PPdBd08E1ie5PSqWjmrpkvpr/ZesmSJkwxK0ogM8+2m1yX5ZpIHgcuA729HQEB3tfbiJEck2ZfuKGT5pher6jtVNa+qFlXVIuBGYLOAkCSNzzDDTe+g+0bS84CvA3+6PTuqqseAc4Dr6M5jXFVVa5JclOT07dmmJGm0hvkK7HeralW//I4k2/3V1KpaQXeye7Dtwkbfk7Z3P5KknWOYkDg0yTK63/6/Cuwz2pIkSdNimJD4E+AY4GzgBcBTk6ygu37ijqq6YoT1SZImaJhvNz3hnhFJFtCFxguA0wBDQpJ2U9s0CyxAVc3QfZ11xdb6SpJ2bc4CK0lqMiQkSU2GhCSpyZCQJDUZEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLUZEhIkpoMCUlSkyEhSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQkSU2GhCSpyZCQJDUZEpKkprGGRJKlSdYlWZ/k/DlePzfJ2iR3JLk+ybPGWZ8k6YnGFhJJ9gIuAU4FjgbOSnL0rG6rgCVVdQxwNfDucdUnSdrcOI8kTgDWV9XdVfUIcCVwxmCHqrqhqn7Yr94ILBhjfZKkWcYZEocB9w2sz/RtLW8EPj3XC0mWJVmZZOXGjRt3YomSpEHjDInM0VZzdkxeAywB3jPX61V1aVUtqaol8+fP34klSpIG7T3Gfc0ACwfWFwD3z+6U5BTg7cAvVdWPx1SbJGkO4zySuBlYnOSIJPsCZwLLBzskOR74EHB6VW0YY22SpDmMLSSq6jHgHOA64C7gqqpak+SiJKf33d4DPBX4eJLbkixvbE6SNAbjHG6iqlYAK2a1XTiwfMo465EkbZlXXEuSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLUZEhIkpoMCUlSkyEhSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQkSU2GhCSpyZCQJDUZEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLUNNaQSLI0ybok65OcP8frT07ysf71m5IsGmd9kqQnGltIJNkLuAQ4FTgaOCvJ0bO6vRF4qKqeC7wX+LNx1SdJ2tw4jyROANZX1d1V9QhwJXDGrD5nAH/VL18NvCJJxlijJGnAOEPiMOC+gfWZvm3OPlX1GPAd4BljqU6StJm9x7ivuY4Iajv6kGQZsKxf/X6SdTtY2zzgmzu4jXHarnozmcG73eK9ndB7t7X97hbv7c4ygr+jXen93Z5anzVMp3GGxAywcGB9AXB/o89Mkr2BpwMPzt5QVV0KXLqzCkuysqqW7KztjdquVO+uVCvsWvXuSrWC9Y7SKGsd53DTzcDiJEck2Rc4E1g+q89y4HX98quBv6+qzY4kJEnjMbYjiap6LMk5wHXAXsBHqmpNkouAlVW1HPgw8NEk6+mOIM4cV32SpM2Nc7iJqloBrJjVduHA8o+AXx9nTb2dNnQ1JrtSvbtSrbBr1bsr1QrWO0ojqzWO5kiSWpyWQ5LUtMeHxNamCpmkJB9JsiHJ6oG2g5N8NsnX+uefmWSNg5IsTHJDkruSrEny1r596mpO8pQkX0lye1/rf+jbj+inhPlaP0XMvpOudVCSvZKsSnJtvz619Sa5J8mdSW5LsrJvm7rPAkCSg5JcneSr/ef3pVNc65H9e7rp8d0kbxtVvXt0SAw5VcgkXQ4sndV2PnB9VS0Gru/Xp8VjwHlVdRTwEuAt/fs5jTX/GDi5qo4FjgOWJnkJ3VQw7+1rfYhuqphp8lbgroH1aa/35VV13MDXM6fxswDwX4G/q6rnAcfSvcdTWWtVrevf0+OAFwE/BD7JqOqtqj32AbwUuG5g/QLggknXNavGRcDqgfV1wKH98qHAuknXuIXaPwW8ctprBvYHbgVeTHdB0t5zfT4m/aC7tuh64GTgWrqLT6e53nuAebPapu6zADwN+Cf6c7TTXOsctf8y8KVR1rtHH0kw3FQh0+afVdUDAP3zIROuZ079DL7HAzcxpTX3Qze3ARuAzwL/G/h2dVPCwPR9Ht4H/BHwk379GUx3vQV8Jskt/SwJMJ2fhWcDG4H/1g/lXZbkAKaz1tnOBK7ol0dS754eEkNNA6Jtk+SpwN8Ab6uq7066npaqery6Q/YFdBNQHjVXt/FWNbckvwpsqKpbBpvn6DoV9fZOrKoX0g3nviXJL066oIa9gRcCf1FVxwM/YEqGlrakP/90OvDxUe5nTw+JYaYKmTbfSHIoQP+8YcL1PEGSfegC4q+r6hN981TXXFXfBj5Pdx7loH5KGJiuz8OJwOlJ7qGbQflkuiOLaa2Xqrq/f95AN2Z+AtP5WZgBZqrqpn79arrQmMZaB50K3FpV3+jXR1Lvnh4Sw0wVMm0Gpy55Hd24/1Top3X/MHBXVV088NLU1ZxkfpKD+uX9gFPoTlbeQDclDExJrQBVdUFVLaiqRXSf07+vqrOZ0nqTHJDkwE3LdGPnq5nCz0JV/V/gviRH9k2vANYyhbXOchb/f6gJRlXvpE+8TPoBnAb8I9149NsnXc+s2q4AHgAepftt541049DXA1/rnw+edJ0D9f4C3XDHHcBt/eO0aawZOAZY1de6Griwb3828BVgPd1h/JMnXesctZ8EXDvN9fZ13d4/1mz6tzWNn4W+ruOAlf3n4RrgZ6a11r7e/YFvAU8faBtJvV5xLUlq2tOHmyRJW2BISJKaDAlJUpMhIUlqMiQkSU2GhCSpyZCQJDUZEtpjJHm8n39/dZKPJ9l/B7b1ziR/0C9/eQv9Dkry5u3dz6xt7ZfkC/0U9yQ5PsnFSb44MDWHtFMZEtqTPFzdPPzPBx4B/u3gi+ls87+Jqvr5Lbx8ELBNIbGFOn4b+ERVPd7vd1VVnUt3de2/2ZZ9SMMyJLSn+l/Ac5Ms6u9E9kG6e0osTPKa/q51tyX50MBv7m9PdxfDzwGb5vkhyfcHln8ryR3p7nj3UeBdwHP6bb2n73NufzSzOsnb+rbN6pij5rMZmI8nyf9IchLdNBJn79R3R+p5iKo9Tj80cyrwd33TkcAbqurNSY6i+638xKp6tP9P++wka+gm1jue7t/NrcAts7b7c8Db+5/9ZpKD6W5o8/zqpiQnyYuAN9Dd4CjATUm+QHdXuZ/WMUfN+wLPrqp7BpqPoZtr6DvAv9jBt0WakyGhPcl+/U2GoDuS+DDwz4F7q+rGvv0VdLeEvLmb1Jb96KZcPhj4ZFX9ECDJXLMFnwxcXVXfBKiqB5M8bVafX+i384N+O58AXkY3g+dgHbPNA769aaUPjQOq6sF+/ZEkB1bV94Z7K6ThGBLakzy86Tf6Tfog+MFgE/BXVXXBrH5vY+s39MmQfVp+sIXXHgaeMrB+NE+81/WTgR9tZd/SNvOchPRE1wOvTnIIQJKDkzwL+CLwqv4bRgcC/7Lxs7+R5Bmbfhb4HnDgQJ8vAv8qyf79fRZeRXdUs0VV9RCwV5JNQXEs3TTc9PvbWFWPbvsfV9oyjySkAVW1Nskf092b+Ul09/J4S1XdmORjdPfIuJc5/mOvqjVJ/hPwhSSPA6uq6vVJvpRkNfDpqvrDJJfT3QMC4LKqWpXunuBb8xm64arP0YXEpqGplwMrtvOPLG2R95OQdhFJjgfOrarXJrkVeFVV3duf17igqtZNuETthhxuknYRVbUKuDHJKuCzfUDsC1xjQGhUPJKQJDV5JCFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLUZEhIkpoMCUlS0/8DjG3PrqKhI8kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(np.arange(data.shape[1]), data.mean(axis=0))\n",
    "plt.xlabel(\"Predictor ($j$)\")\n",
    "plt.ylabel(r\"$P(j \\in M \\mid \\mathbf{Y}, \\mathbf{X})$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the **actual predictors** have very high posterior probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([13, 31, 38, 55]),)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predicted_model = data.mean(axis=0) > 0.7\n",
    "predicted_model.nonzero()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uncomment the following line to see number of times different models were sampled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pd.DataFrame(data).assign(count=-1).groupby(list(range(70))).count()"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}