main_test_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test over Synthetic Data - part 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll try this method on another boolean function that has more predictors (both actual and false).\n",
    "\n",
    "We will generate 1500 data points. Each data point will contain 100 predictors ($X_i = x_{i1}, x_{i2}, \\dots, x_{i100}$) and the response of the system ($y_i$). To generate the random data points $\\mathbf{X}=<X_1, X_2, \\dots, X_{1500}>$, we sample a thousand times from the binomial distribution with $n=100$ and $p=0.5$."
   ]
  },
  {
   "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",
    "    sys.path.append('noisy-generalized-boolean')\n",
    "    project_root = \"noisy-generalized-boolean/\""
   ]
  },
  {
   "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\n",
    "import os"
   ]
  },
  {
   "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 = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 1., 0., ..., 1., 1., 1.],\n",
       "       [1., 1., 1., ..., 0., 0., 1.],\n",
       "       [1., 0., 1., ..., 0., 0., 1.],\n",
       "       ...,\n",
       "       [0., 0., 1., ..., 1., 1., 1.],\n",
       "       [0., 0., 1., ..., 0., 0., 1.],\n",
       "       [0., 1., 1., ..., 0., 0., 0.]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.random.rand(1500, 100)\n",
    "X = np.round(X)\n",
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have $\\mathbf{X}$, it's time to generate our synthetic system's response $\\mathbf{Y}$. We say that the underlying boolean function depends only on the following 6 predictors out of 100 predictors:\n",
    "                $$M = <x_{13}, x_{31}, x_{38}, x_{55}, x_{77}, x_{89}>$$\n",
    "\n",
    "We'll represent this function by $F: \\{0,1\\}^6 \\rightarrow \\mathbb{Z}_{4} $, where $\\mathbb{Z}_{4} = \\{0, 1, 2, 3\\}$.\n",
    "\n",
    "The choice of $F$ is very simple. If the *boolean OR* of the actual predictors is *True*, then the bool_func generates:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "     y &= 0 \\:\\text{ w.p. }\\: 0.9 \\\\\n",
    "       &= 1 \\:\\text{ w.p. }\\: 0.05 \\\\\n",
    "       &= 2 \\:\\text{ w.p. }\\: 0 \\\\\n",
    "       &= 3 \\:\\text{ w.p. }\\: 0.05\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "And if the *boolean OR* evaluates to be *False*, then bool_func generates\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "     y &= 0 \\:\\text{ w.p. }\\: 0 \\\\\n",
    "       &= 1 \\:\\text{ w.p. }\\: 0.1 \\\\\n",
    "       &= 2 \\:\\text{ w.p. }\\: 0.7 \\\\\n",
    "       &= 3 \\:\\text{ w.p. }\\: 0.2\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Here \"w.p.\" abbreviation stands for *with probability*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bool_func(xi):\n",
    "    \"\"\"Takes an row of truth table and returns the response \n",
    "    according to the boolean function\"\"\"\n",
    "    on_true = [0.9, 0.05, 0, 0.05]\n",
    "    on_false = [0, 0.1, 0.7, 0.2]\n",
    "    weights =  on_true if np.any(xi) else on_false\n",
    "    return np.random.choice(4, p=weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the above defined function, we'll generate the response of the system for input $\\mathbf{X}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1500, 6)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_actual = X[:, [89, 77, 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 $1500$ rows in $\\mathbf{X}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1500,)"
      ]
     },
     "execution_count": 7,
     "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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = []\n",
    "if not load_samples:\n",
    "    samples = main.metropolis(X, Y, q=4, iterations=50000, interval=5000, show_progress=True)\n",
    "    pd.DataFrame(samples).to_csv(f\"{project_root}data/mcmc-2.csv\")\n",
    "else:\n",
    "    samples = pd.read_csv(f\"{project_root}data/mcmc-2.csv\", index_col=0).values.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a href='data/mcmc-2.csv' target='_blank'>data/mcmc-2.csv</a><br>"
      ],
      "text/plain": [
       "/home/rohan/Projects/PythonProjects/generalized-noisy-boolean/data/mcmc-2.csv"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FileLink(f\"{project_root}data/mcmc-2.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Click on the link above to download the generated samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ioff() # turn off interactive mode for matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(49581, 100)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = np.array(samples)\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAH/RJREFUeJzt3XucHGWd7/HPb3puuYckkwQTMEQiICgQRxSBcxBQo6iwrLJeNyorZ3XX+zka1/uungXd46qLLOQgElbkIhdhkYssgoBAwiQkJJBAQggkJCRDMrlnMrff/lFPT7on3dPVk+nLTH3fr9e8uqq6quv39FT3r5/nqXrK3B0REZG0mkoHICIi1UWJQUREsigxiIhIFiUGERHJosQgIiJZlBhERCSLEoOIiGRRYhARkSxKDCIikqW20gEMxKRJk3zGjBmVDkNEZEhZvHjxq+7eVGi9IZkYZsyYQUtLS6XDEBEZUszsxTjrqSlJRESyKDGIiEgWJQYREcmixCAiIlmUGEREJEtZz0oys3XALqAb6HL3ZjObANwIzADWARe6e1s54xIRkQMqUWN4h7uf5O7NYX4ecL+7zwLuD/MiIlIh1dCUdB6wIEwvAM4v1Y5uXbKB6xbGOo1XRCSxyp0YHPiDmS02s4vDsinuvgkgPE7OtaGZXWxmLWbW0traOqCd37FsIzc+sX5A24qIJEW5r3w+zd03mtlk4D4zWxV3Q3efD8wHaG5u9oHs3AAf0JYiIslR1hqDu28Mj1uA24BTgM1mdjhAeNxSqv2bGY4yg4hIf8qWGMxslJmNSU8D7wJWAHcAc8Nqc4HbSxVDjanGICJSSDmbkqYAt5lZer+/cfd7zOwJ4CYzuwh4CfhQ6UIwepQYRET6VbbE4O5rgRNzLN8KnF2OGMzAVWUQEelXNZyuWjZW6QBERIaAZCUG9TGIiBSUrMSAzkoSESkkWYlBNQYRkYISlRhqzFRfEBEpIFGJAYMeVRlERPqVqMRggKoMIiL9S1ZiUFOSiEhByUoM6AI3EZFCkpUYTC1JIiKFJCsxoNNVRUQKSVZi0LDbIiIFJSwxqMYgIlJIshIDpsQgIlJAshKDht0WESkoWYkBnZUkIlJIshKD+hhERApKVmLQsNsiIgUlKzGoxiAiUlDCEoPGShIRKSRhiUFnJYmIFJKsxICakkRECklWYtAgeiIiBSUrMWBqShIRKSBZiUE1BhGRghKVGDq7ne17OysdhohIVUtUYrh+0UuVDkFEpOolKjGIiEhhSgwiIpJFiUFERLIoMYiISJayJwYzS5nZk2Z2Z5g/yswWmtlqM7vRzOpLte+LTj+qVC8tIjJsVKLG8EVgZcb8pcC/uvssoA24qFQ7HlWfKtVLi4gMG2VNDGY2HTgXuCrMG3AWcHNYZQFwfgkDKNlLi4gMF+WuMfwU+BrQE+YnAtvdvSvMbwCmlTkmERHJULbEYGbvA7a4++LMxTlWzTlqhZldbGYtZtbS2to6sBgGtJWISLKUs8ZwGvABM1sH3EDUhPRTYLyZ1YZ1pgMbc23s7vPdvdndm5uamg4pEA2kJyKSX9kSg7t/w92nu/sM4MPAH939Y8ADwAfDanOB20sVg7oYREQKq4brGL4OfMXM1hD1Ofyy1DtUhUFEJL/awqsMPnd/EHgwTK8FTinHfk29DCIiBVVDjaHsVGEQEckvUYlBfQwiIoUlKjH85L7nAHhg1ZYKRyIiUr0SlRjSbmpZX+kQRESqViITg/oYRETyS2Zi0PmqIiJ5JTIx9CgviIjklcjEoJOTRETyS2RiEBGR/JQYREQkixKDiIhkUWIQEZEsiUwMOilJRCS/RCYGERHJr+jEYGajzCxVimDKZfzIukqHICJStQomBjOrMbOPmtnvzWwLsArYZGZPm9mPzWxW6cMcHJd/bDYA73/TayociYhI9YpTY3gAeB3wDWCqux/h7pOBM4DHgUvM7OMljHHQvGb8iEqHICJS9eLcwe0cd+/su9DdtwG3ALeY2ZBqm3F1P4uI5FWwxpArKQxknWqwZ38XAI+s3lrhSEREqlfBGoOZHRnztba7+85DjKekVm/eBcDVf36B77z/DRWORkSkOsVpSloQYx0HrgGuPaRoSsx0b08RkYLiJIZ/Apa7e2upgyk15QURkcLiJIb/AraYWQ+wAngKWB4en3H3/SWMb1ApL4iIFBYnMXwB+DRwE/AocAzwZuCTwHHA1FIFN+hUZRARKSjOWUmXAacR9SP8FOgEvuju73D3oZMUUI1BRCSOWENiuPs+d78UOBM4GlhkZm8tZWCloAqDiEhhcU5XPYOoyejY8DgZ2AVMLG1opdXe2U1j3ZAe8klEpCTi9DH8CVgGXA/83N3XlTSiEqrJqDJsaNvL0ZPHVDAaEZHqFCcxfBZ4I3Au8FUz20p0VtJyYIW7/66E8Q2qzJYk16gYIiI5FUwM7n5l5ryZTQfeRJQs/hIYMokhk/KCiEhucWoMWdx9A7ABuGvwwymtzM7nrbs7YErlYhERqVZlGyvJzBqBh4CGsN+b3f27ZnYUcAMwAVgCfMLdO2LusyiZQ2L818rNnPq6Id1/LiJSEnHHSnL6vwwgzlhJ+4Gz3H13GKb7ETO7G/gK8K/ufoOZXQFcBPx7nOCLZXmmRUTkgDh9DO8YjB25uwO7w2xd+HPgLOCjYfkC4HuUKDGIiEhhcW7tWfDHdZx1wnopM1sKbAHuA54naoLqCqtsAKbFea2BqEsdKK4udhMRyS3WrT3N7PN9+xrMrN7MzjKzBcDcODtz9253PwmYDpxCdMHcQavl2tbMLjazFjNraW0d2ECvZx83uXe6pkaZQUQklziJYQ7QDVxvZhvN7BkzWwusBj5C1D9wTTE7dfftwIPA24DxZpZu0poObMyzzXx3b3b35qampmJ216u25kBxU6oyiIjkFKePoR24HLg8dBpPAvaFL/fYzKwJ6HT37WY2AjgHuBR4APgg0ZlJc4HbiytCMTEcmK5VjUFEJKeirmMI93beNMB9HQ4sMLMUUU3lJne/08yeAW4wsx8ATwK/HODrFyVVE2v8QBGRxCn6AreBcvengJNzLF9L1N9Qcpk1hpOPHF+OXYqIDDmJ+tlsGVcv1KiPQUQkp9iJwSIfN7PvhPkjzawsv/QHi3KBiEhhxdQYLgdOJToTCaJ7Mvxi0CMqIeUFEZHCiuljeKu7zzazJwHcvc3M6ksUV0lkXofXo3G3RURyKqbG0BnOKHLoPf20pyRRlYHSgohIbsUkhp8DtwGTzeyHwCPAP5ckqhLJvlGPUoOISC6xm5Lc/TozWwycTfQde767ryxZZCWQ2fmsvCAiklsxZyXdD8x091+4+2XuvtLM5pcwtkGX2cfgakwSEcmpmKako4Cvm9l3M5Y1D3I8ZdMzZHtHRERKq5jEsJ2oGWmKmf2nmY0rUUxlofqCiEhuxSQGc/cud/8ccAtR5/PkAttUrbY9Jbl7qIjIkFfMdQxXpCfc/RozWw783eCHVB4j6lOVDkFEpCoVc1bSlX3mFwOfHvSIykRNSSIiucW5tecj4XGXme0Mf7vS86UPsTR0HYOISG5xbtRzengcU/pwykdDYoiI5FbMdQwfMrMxYfpbZnarmR10f4WhQqeriojkVsxZSd92911mdjrwbmABGR3SQ43qCyIiuRWTGLrD47nAv7v77cCQGl01k5qSRERyKyYxvGxmVwIXAneZWUOR21cVdT6LiORWzBf7hcC9wBx33w5MAP5PSaIqgx7lBRGRnIq5jmEvcGvG/CZgUymCKgdVGEREchuyTUGHSn0MIiK5JTYxqI9BRCS3xCWGT759BqA+BhGRfOIMiZFzKIyhOiTGF8+eBagpSUQknzhDYgyroTBqwl3clBdERHIrZkgMM7OPm9m3w/wRZnZK6UIrDQslVo1BRCS3YvoYLgdOBT4a5ncDvxj0iEpMNQYRkf4Vc6Oet7r7bDN7EsDd28xsyA2JURPlBdUYRETyKKbG0GlmKcL4c2bWBAy5MUqNKDPorCQRkdyKSQw/B24DJpvZD4nu+fx/SxJVCYWWJLbu3l/ZQEREqlQxQ2JcZ2aLgbMBA85395Uli6xE0n0MoxqKaUUTEUmOor4d3X0VsGogOzKzI4BrgalETVDz3f1nZjYBuBGYAawDLnT3toHsI45U6GRIP4qISLaCicHMvtLf8+7+k5j76gK+6u5Lwp3gFpvZfcAngfvd/RIzmwfMA74e8zWLls4H3epkEBHJKU6NIX2B2zHAW4A7wvz7gYfi7ihzNNZwJ7iVwDTgPODMsNoC4EFKmBis93RVJQYRkVziXPn8fQAz+wMw2913hfnvAb8dyE7NbAZwMrAQmBKSBu6+ycwm59nmYuBigCOPPHIgu+2VqjGdlSQikkcxZyUdCXRkzHcQ9QsUxcxGA7cAX3L32GMtuft8d2929+ampqZid5ulxnQdg4hIPsV0Pv8HsMjMbgvz5xM1/cRmZnVESeE6d0/f9GezmR0eaguHA1uKec2BMFONQUQkn9g1Bnf/IfApoA3YBnzK3f857vYWNe7/EljZp8P6DmBumJ4L3B73NQeqxtTHICKST7En83cTnWrqFH/V82nAJ4DlZrY0LPsH4BLgJjO7CHgJ+FCRr1u0GjOdlSQikkfsxGBmXwQ+Q9QUZMCvzWy+u/9bnO3d/ZGwXS5nx41jMKTUlCQiklcxNYaLiAbS2wNgZpcCjwGxEkM1MXU+i4jkVcxZSUbUlJTWTf4aQFWrqTH1MYiI5FFMjeFXwMJwVpIRXZh2dUmiKrEaNSWJiORVzCB6PzGzB4HTw6JPuvvSfjapWrqOQUQkvzhjJd3Rd1F4PMfMcPcPDH5YpRVdx6DEICKSS5waw6nAeuB6oiEshmS/QiYDtu7uKLieiEgSxUkMU4F3Ah8hut/z74Hr3f3pUgZWSlv3dFCXKqbfXUQkOQp+O7p7t7vf4+5zgbcBa4AHzezzJY+uRKaObaS+VolBRCSXWJ3PZtYAnEtUa5hBdJvPW/vbppo11NbQpdOSRERyitP5vAA4Abgb+L67ryh5VCWWqjG6e4od0UNEJBni1Bg+AewBXg98IX2jG6I+XHf3sSWKrWRSNUZXt2oMIiK5xLlRz7BrjK9NaRA9EZF8ht2XfhypGvUxiIjkk8jEUFtjtO3VdQwiIrkkMjG07elQU5KISB7F3qhnWBg/so7amkTmRBGRghL57dhYl8JRjUFEJJdEJgYNuy0ikl8iE4Pu4CYikl8iE4NqDCIi+SU0MaBbe4qI5JHIxKAb9YiI5JfIxBDVGCodhYhIdUpkYjD1MYiI5JXIxKA+BhGR/BKaGNTHICKST4ITQ6WjEBGpTolMDOgCNxGRvBKZGGrM0FBJIiK5JTQxqMYgIpJPQhOD+hhERPIpW2Iws6vNbIuZrchYNsHM7jOz1eHxsPLEohqDiEg+5awxXAPM6bNsHnC/u88C7g/zJVdjpiufRUTyKFticPeHgG19Fp8HLAjTC4DzyxGLAfu7usuxKxGRIafSfQxT3H0TQHicXI6dtnf18OrujnLsSkRkyKl0YojNzC42sxYza2ltbT2k1xpZl6LGBikwEZFhptKJYbOZHQ4QHrfkW9Hd57t7s7s3NzU1HdJOJ46uJ6XMICKSU6UTwx3A3DA9F7i9HDvV6aoiIvmV83TV64HHgGPMbIOZXQRcArzTzFYD7wzzJacL3ERE8qst147c/SN5njq7XDH00umqIiJ5VbopqSLS3Qu6J4OIyMESmhiizKB+BhGRgyU0MUSP6mcQETlYIhOD9dYYlBhERPpKZGJINyUpL4iIHCyRicHUlCQiklciE8OBs5IqG4eISDVKaGJQH4OISD6JTAym01VFRPJKZGLQBW4iIvklNDGoxiAikk8iE4POShIRyS+hiUHXMYiI5JPIxKA+BhGR/BKZGNJ9C/u7eiobiIhIFUpkYugICUEVBhGRgyUyMUwYVQdAtzKDiMhBEpkYdOWziEh+yU4MupBBROQgiUwMqXBakpqSREQOlsjEcKDGUOFARESqUEITQ/SoPgYRkYMlMjH0NiWpj0FE5CCJTAw16mMQEckrkYkh1TtWkhKDiEhfiUwM6c7nbnU+i4gcpLbSAVRCTUiHF175GE99710s37CDj121EIBvnXscf3PGzKz1N27fx9sv+SMAd37+dN73b49QX1tDR1cP9akaxo2so3XX/oP2c8asSWzf20nb3g42tO0D4JQZE1i0bhsAI+tT/PGrZzJ1XCOfubaF+57ZDMC4EXUcPXk0i19sA6AuZUwc1cArO9u58/Onc8K0cXR19/CxqxZy5jGT+eyZrxv8N6mPlZt28p6fPczYxloe/tpZjBtZl3fdXz7yAv905zMHLa9P1fCD80/gwrcckXfbb9y6nBqDH/7FG/uN54Fnt/CpXz3ROz+2sZZ/PO8EvnTj0rzbfPStR3LL4g15x8j6+pxj+d2TL/Ps5l1Zy9ddcm6/sQyG9s5uPvL/H+eN08bxj+ed0O+6G7fv49PXPMFnzpjJmMZaLv6PxQAcPXk0c46fSsuL2/j1RW/le//5NNv2dHDX8lcK7v9Tp83gu+8/Pu/zM+b9HoDXjGvkgtnTueyBNbHLdvTk0Vz5iTdz1cMvcP2il/Ku9/8+dCJ/+ebpBV+vu8e58MrHWPxiGz/78Emcd9K0gtvs3t/FBZf/mX/50Im8afr4guv/6blWvnv7Cua951jmnHA4X7lxKd3u3L50IwDHTh3DOcdN4bIH1jC2sZad7V0ANNbVsPQ77+Kbt61g4Qtb+da5xzHnhMMBuGv5Jj533RKu+PhsfnjXStZvi74TLpg9jZ9ceFLBmMopkTUGw3qn12/by9dveap3/ge/X3nQ+jcv3tA7fd3C6MBOj7fU0d2TMykAPLz6VZa/vKM3KQC9SQFgb0c3a7bsBuhNCgA79nX2JgWAzm7nlZ3tQPSlm15n4QvbuPSeVYWKOyhuCe/BzvYuXti6p991cyUFiN6rr2W817lcv+il3ve4Pz/os4+d7V39JgWA3yx8qd+BEy+9Z9VBSaFctuzcz5Mvbefax14suO5zm3ex6pVdXPGn5/lyRpnXbNnNZQ+s4fG129jZ3sWvH38pVlIA+NWf1+V9LvNC0I072otKCum4Vm3a1W9SAPjqb5fFer2dGZ+Pb/9uRaxtlrzYxnObd/Pje5+NtX7Lum2s27qXRS9E+7n1yZd7kwLAqld29b4P6aQA0N4ZfR/84ZlX2NC2j4UvHPi8f+66JQD87a+X9CYFgFuXvBwrpnJKZGIYN+LAr92Orh7qUv2/DZmntXYM8oisHd3dRZ0dta+jO2xXuXawajibq71zeLUDdnR3x183HIM97uzpyL3dYB6nnYNwwU8x5Sv8Wgfi2dcZ73XTn+H0vVgK7qP3h1/xce/v6un9jAz290W5JDIxNNQdKPb+rp7e01fjGOwv5I6unqIOnvT+91fwi3GoHuzVrJhEF+cYHMz/UVf3of8QGMzjNbNscc8fSa8X96OerlkO5H3c39V9SNtXAxuKZ+Y0Nzd7S0vLgLdfv20vZ/zoAQCmjR/By9v3ZT0/a/LorPnVobkHogNrMH8wTx3byMiGFGtb+2+eyTRr8mg6unt4ceve3vlSy3wPpo0fwcj6VKx1c+kv3vS2hcpUaB+DqRzvb3tXd2/zQqH97WzvZPPO3M2Xaa+dOLL3+Igr33673Ys6PnNpGtOQt8k1TgyZMo/9uNusb9vbm3zjrL95Zzs727sY01DL1HGNRR1vR0wY0fu/TG8P/R+zxRxjv5z7Fo6cODL2+pnMbLG7Nxdaryo6n81sDvAzIAVc5e6XlHJ/08aPYNLoel7d3cGJR4xjytgGlry0HYg612ZNyf4nzWwaxb1PR30A7z5+KneveIW6lNHZ7ZwwbSwrXt7Zu+70w0b09imkk87I+hR7Q5X/DYePZd3WPYwfUcfGHe3Mfm3UEXb8a8axdH0b67ftI1VjjB9Rx9Y9Hb2v2/zaw2h5sY1zjptCfW30s2dXexdHN41m0pj6Er1TBxw9eTR3r4jaq088Yly/63b3OGtfzf1FcuzUMcxsGpV3287uHszsoP9BX1PGNvLImlezlp185HieDP/HvprGNDBz0qisNt/XTxnN3LfP4PdPbeLR57dyxqxJPLz61awTBICCsQyWEXUppoxtZExj4Y/lw8+9ytuPnsgrO9pZtmEHr2saxRETRrJnfxdPrGvj+NeM5dVd+0nVWFYb+GEj66ivrWFfRzfHTh1Lbcp49PmtTBnb0G85MxPDzEmj8v5/+7pg9jQaalPs2NfBsvU7en+EpT8/jXU1vV/Y//P1TYxqyP+DI9PYxjqWv7yDU2dO5LBR+U+ESHtd02juefoV3jZzAhNGFf68zJoymuc27+b14T1ZvWU3k0bXM2l0Ax9883TWbd3DvU9vpnXX/t4TUUbWp+jo6uGN08Zx4vTxTBxVT+vuA8lwZEMty9Zv5x3HNHH05NFc8+g6Orud+lRNUcdYfW3pG3oqXmMwsxTwHPBOYAPwBPARd8/dg8mh1xhERJIobo2hGvoYTgHWuPtad+8AbgDOq3BMIiKJVQ2JYRqwPmN+Q1gmIiIVUA2JIdd5Age1b5nZxWbWYmYtra2tZQhLRCSZqiExbAAyL4WdDmzsu5K7z3f3ZndvbmpqKltwIiJJUw2J4QlglpkdZWb1wIeBOyock4hIYlX8dFV37zKzvwfuJTpd9Wp3f7rCYYmIJFbFEwOAu98F3FXpOEREpDqakkREpIpU/AK3gTCzVqDwMJS5TQJeLbjW8JLEMkMyy53EMkMyyz2QMr/W3QuevTMkE8OhMLOWOFf+DSdJLDMks9xJLDMks9ylLLOakkREJIsSg4iIZEliYphf6QAqIIllhmSWO4llhmSWu2RlTlwfg4iI9C+JNQYREelHohKDmc0xs2fNbI2Zzat0PMUys6vNbIuZrchYNsHM7jOz1eHxsLDczOznoaxPmdnsjG3mhvVXm9ncjOVvNrPlYZufW9wb5JaQmR1hZg+Y2Uoze9rMvhiWD9tym1mjmS0ys2WhzN8Py48ys4Uh/hvDEDKYWUOYXxOen5HxWt8Iy581s3dnLK/Kz4KZpczsSTO7M8wnoczrwvG31MxawrLKHt/unog/ouE2ngdmAvXAMuANlY6ryDL8D2A2sCJj2Y+AeWF6HnBpmH4vcDfR6LVvAxaG5ROAteHxsDB9WHhuEXBq2OZu4D1VUObDgdlhegzRTZ3eMJzLHeIYHabrgIWhLDcBHw7LrwA+G6Y/B1wRpj8M3Bim3xCO8wbgqHD8p6r5swB8BfgNcGeYT0KZ1wGT+iyr6PGdpBrDkL8hkLs/BGzrs/g8YEGYXgCcn7H8Wo88Dow3s8OBdwP3ufs2d28D7gPmhOfGuvtjHh1N12a8VsW4+yZ3XxKmdwErie7XMWzLHWJP3yC4Lvw5cBZwc1jet8zp9+Jm4Ozwq/A84AZ33+/uLwBriD4HVflZMLPpwLnAVWHeGOZl7kdFj+8kJYbhekOgKe6+CaIvUWByWJ6vvP0t35BjedUIzQUnE/2CHtblDk0qS4EtRB/y54Ht7p6+gXNmnL1lC8/vACZS/HtRaT8Fvgb0hPmJDP8yQ5T0/2Bmi83s4rCsosd3VQyiVyaxbgg0jOQrb7HLq4KZjQZuAb7k7jv7aSYdFuV2927gJDMbD9wGHJdrtfBYbNly/SCsaJnN7H3AFndfbGZnphfnWHXYlDnDae6+0cwmA/eZ2ap+1i3L8Z2kGkOsGwINQZtDdZHwuCUsz1fe/pZPz7G84sysjigpXOfut4bFw77cAO6+HXiQqD15vJmlf8xlxtlbtvD8OKImx2Lfi0o6DfiAma0jauY5i6gGMZzLDIC7bwyPW4h+BJxCpY/vSne8lOuPqHa0lqhDKt35dHyl4xpAOWaQ3fn8Y7I7qX4Ups8lu5NqkR/opHqBqIPqsDA9ITz3RFg33Un13ioorxG1i/60z/JhW26gCRgfpkcADwPvA35Ldkfs58L035HdEXtTmD6e7I7YtUSdsFX9WQDO5EDn87AuMzAKGJMx/Sgwp9LHd8UPgjL/E95LdFbL88A3Kx3PAOK/HtgEdBL9EriIqF31fmB1eEwfDAb8IpR1OdCc8TqfJuqUWwN8KmN5M7AibHMZ4QLICpf5dKKq71PA0vD33uFcbuBNwJOhzCuA74TlM4nOMFkTvjAbwvLGML8mPD8z47W+Gcr1LBlno1TzZ4HsxDCsyxzKtyz8PZ2Oq9LHt658FhGRLEnqYxARkRiUGEREJIsSg4iIZFFiEBGRLEoMIiKSRYlBEsfMdofHGWb20UF+7X/oM//oYL6+SDkoMUiSzQCKSgxmliqwSlZicPe3FxmTSMUpMUiSXQKcEcbB/3IYuO7HZvZEGOv+fwGY2ZkW3RPiN0QXFWFmvwuDnj2dHvjMzC4BRoTXuy4sS9dOLLz2ijA2/l9lvPaDZnazma0ys+vS4+Wb2SVm9kyI5V/K/u5IYiVpED2RvuYB/9vd3wcQvuB3uPtbzKwB+LOZ/SGsewpwgkdDOQN82t23mdkI4Akzu8Xd55nZ37v7STn2dQFwEnAiMCls81B47mSioRw2An8GTjOzZ4C/AI51dw+D6YmUhWoMIge8C/jrMNz1QqJhCWaF5xZlJAWAL5jZMuBxosHLZtG/04Hr3b3b3TcDfwLekvHaG9y9h2jIjxnATqAduMrMLgD2HnLpRGJSYhA5wIDPu/tJ4e8od0/XGPb0rhQNC30OcKq7n0g0rlFjjNfOZ3/GdDdQ69E9Bk4hGlX2fOCeokoicgiUGCTJdhHdLjTtXuCzYZhvzOz1ZjYqx3bjgDZ332tmxxKNXJnWmd6+j4eAvwr9GE1Et2ldlC+wcP+Jce5+F/AlomYokbJQH4Mk2VNAV2gSugb4GVEzzpLQAdxK7tsg3gP8rZk9RTSC5+MZz80HnjKzJe7+sYzltxHdd3cZ0WixX3P3V0JiyWUMcLuZNRLVNr48sCKKFE+jq4qISBY1JYmISBYlBhERyaLEICIiWZQYREQkixKDiIhkUWIQEZEsSgwiIpJFiUFERLL8N4aSmAoRJ2rqAAAAAElFTkSuQmCC\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 10000 iterations. We discard the samples before that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(39581, 100)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = data[10000:]\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next plot gives the posterior probability that predictor $j$ is a part of the model given $Y$ and $X$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "     \\begin{align*}\n",
    "     P(j \\in M \\mid \\mathbf{Y}, \\mathbf{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": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEOCAYAAABxdpuaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFyVJREFUeJzt3X20XXV95/H3x2B8QCtPcWoJMVAzCoqCRrRarFVUFAs64oijFquWVZWpDjqdWKfY0mUL1bHWJbawkKlju0BFa7NKLCKido0DJhAQAlJjysOdMBIFH0Eg+J0/9g6eXO4l55fcc+7Nve/XWnfd/fDb53zP3cn93P3be/92qgpJklo8ZLYLkCTtfgwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnN9pjtAkZlv/32q+XLl892GZK0W7niiiu+V1VLdtRu3obH8uXLWbdu3WyXIUm7lSQ3DdPObitJUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDC8LyVReyfNWFs12GNG8YHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySp2R7jfLMkRwN/BSwCzqmq0yetPwV4C7AV2AK8qapu6tfdB1zTN725qo4dW+HSHDL4ON0bTz9mFivRQja28EiyCDgTeBEwAaxNsrqqrhtoth5YWVV3Jnkr8BfAa/p1d1XVYeOqV5I0vXF2Wx0BbKyqTVV1D3A+cNxgg6q6tKru7GcvA5aOsT5J0pDGGR77A7cMzE/0y6bzZuALA/MPT7IuyWVJXjHVBklO6tus27Jly65XLEma0jjPeWSKZTVlw+T1wErgNwYWL6uqzUkOAr6c5Jqq+s52L1Z1NnA2wMqVK6d8bUnSrhvnkccEcMDA/FJg8+RGSY4C3gscW1V3b1teVZv775uArwCHj7JYSdL0xhkea4EVSQ5Mshg4AVg92CDJ4cBZdMFx28DyvZM8rJ/eD3guMHiiXZI0RmPrtqqqrUlOBi6iu1T33KrakOQ0YF1VrQY+ADwK+EwS+MUluQcDZyX5OV3gnT7pKi1J0hiN9T6PqloDrJm07NSB6aOm2e7rwKGjrU6SNCzvMJckNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNdtjtgtQu+WrLrx/+sbTj5nFSiQtVB55SJKaGR6SpGZjDY8kRye5IcnGJKumWH9KkuuSfDPJJUkeP7DuxCTf7r9OHGfdGp/lqy68/0vS3DW28EiyCDgTeClwCPDaJIdMarYeWFlVTwUuAP6i33Yf4H3As4AjgPcl2XtctUuStjfOI48jgI1Vtamq7gHOB44bbFBVl1bVnf3sZcDSfvolwMVVdXtV3QFcDBw9prolSZOMMzz2B24ZmJ/ol03nzcAXWrZNclKSdUnWbdmyZRfLlSRNZ5zhkSmW1ZQNk9cDK4EPtGxbVWdX1cqqWrlkyZKdLlSS9ODGGR4TwAED80uBzZMbJTkKeC9wbFXd3bKtJGk8xhkea4EVSQ5Mshg4AVg92CDJ4cBZdMFx28Cqi4AXJ9m7P1H+4n6ZJGkWjO0O86ramuRkul/6i4Bzq2pDktOAdVW1mq6b6lHAZ5IA3FxVx1bV7Un+lC6AAE6rqtvHVbskaXtjHZ6kqtYAayYtO3Vg+qgH2fZc4NzRVSdJGpZ3mEuSmjWHR5I9+xv+JEkL1A7DI8lDkvynJBcmuQ34FnBrkg1JPpBkxejLlCTNJcMceVwK/CrwHuCXq+qAqnoscCTdXeCn9/dlSJIWiGFOmP9uVf3r5IX9FVA3VdWrkjx0BLVJkuaoYY48rkzy3wbPcyR5dJKPAv8HoKruHVWBkqS5Z5jwmAD+HLgqyZFJXk133uNtwNdHWZwkaW4aJjwOBf6Y7rzHV+hGw10MvKWqfmNklUmS5qwdhkffJfWPwEa6AQpD99yNi0dbmiRprhrmUt0P0g0L8kTgg8DfA0cB1yV512jLkyTNRcN0W50CfAN4RlX9QVW9AXgRcCv9k/4kSQvLMJfqnlRV5wwuqKpLkjwF+MPRlCVJmst2GB6Tg2Ng+T10J9IlSQuMAyNKkprtdHgkeVySh81kMZKk3cOuHHl8EvhWfzWWJGkB2emHQVXVUeke93fIDNYjSdoNDHOfxyOmW1edDTNbkiRprhum2+pbSf7DyCuRJO02hgmPXwY+k+QLSZ4w6oIkSXPfMOc8ngx8CHg5cE2SvwFu2bayqj40otokSXPUMDcJbgSOTXIm8FbgHYOr6YJFkgQsX3Xh/dM3nn7MLFYyWjsMjySHAR+me+zsj+kGR7x5xHVJkuawYbqt1tGdG/kUcEpV3TrakiRJc90w4XEDcHJVXTrqYiRJu4dhwuNpVbV15JVIknYbwzxJ0OCQJG1nrKPqJjk6yQ1JNiZZNcX65yW5MsnWJMdPWndfkqv6r9Xjq1qSNNlOj23VKski4Ey6pxBOAGuTrK6q6waa3Qy8EXj3FC9xV1UdNvJCJUk7NMyluucBj3+wJnTDXD1nBy91BLCxqjb1r3s+cBxwf3hU1Y39up/vqC5J0uwZ5ibB187Qe+3PwJ3pdEcfz2rY/uFJ1gFbgdOr6vMzVJckqdHYuq3ojlAmq4btl1XV5iQHAV9Ock1VfWe7N0hOAk4CWLZs2c5XKkl6UOM8YT4BHDAwvxTYPOzGVbW5/74J+Apw+BRtzq6qlVW1csmSJbtWrSRpWsM8z+MFSWbiN/FaYEWSA5MsBk4AhrpqKsne2x55m2Q/4LkMnCuRJI3XMN1WXwJu609iXwt8E7im/35dVd09zBtV1dYkJwMXAYuAc6tqQ5LTgHVVtTrJM4F/APYGfivJn1TVk4GDgbP6Gh5Cd87D8JCkWTJMePw+8Cbg08DXgScCz6C7pPZguud9DKWq1gBrJi07dWB6LV131uTtvg4cOuz7SJJGa5g7zD9K101UdKPr3gu8o6p+s6qGDg5J0vwx1Anzqrqrqs4Ang88AfhGkpbLbCVJ88gwNwkeSdc99aT++2Ppnuux72hLkyTNVcOc8/gqcDVwHvCRbXeBS5IWrmHC4610J6uPAd6V5Pt0V1tdA1zrnd6StPAMMzzJWYPzSZYCT6ULlFcBhockLTDNw5NU1QTd3eJrdtRWkjQ/jfV5HpKk+cHwkCQ1Gyo8knx+YPqfplouSVo4hj3y2H+IaUnSAjFseNQQ05KkBcJzHpKkZkOHR5JHbpvs5x81kookSXPesOFxNvDRJAGq//6RfrkkaYEZdlTdc+geBLUWOAi4ArimXy5JWmCGvsO8qj6U5MPAEmBLVf18dGVJkuaypuFJ+sD47ohqkSTtJrzaSpLUbIfhkeTP+u/PS7LX6EuSJM11w3Rbfbr//rvAU/tLdq+ne0DU1cDlVXXLiOqTJM1Bw4TH1QBV9QaAJA8FDgGeBvwasCrJl6pq1ciqlCTNKcOc87g0yX9Osgygqu6tqquB8+me6XE98JYR1ihJmmOGCY+jgfuA85LcmuS6JJuAbwMnAP8DWDbCGiVJc8wwj6H9GfAx4GN9l9V+wF1V9YNRFydJmpuGudrqxCTfS3I7cA7wE4NDkha2Ybqt/gh4EfAk4Gbgz0ZakSRpzhvmaqsfVdX6fvqPklw+yoIkSXPfMEcej0tyUpIjkywBHrqzb5bk6CQ3JNmY5AGX9vY3Il6ZZGuS4yetOzHJt/uvE3e2BknSrhvmyON9wFOB1wGHAo9Ksobu/o9vVtV5w7xRkkXAmXRdYBPA2iSrq+q6gWY3A28E3j1p2336OlbSPb3win7bO4Z5b0nSzBrmaqvtntmRZCldmBwKvAwYKjyAI4CNVbWpf53zgeOA+8Ojqm7s100esfclwMVVdXu//mK6S4iHfW9J0gxqGlUXoKom6I4c1jRuuj8wOIzJBPCsXdh2/8b3lyTNkHGOqpspltVMbtufm1mXZN2WLVuaipMkDW+c4TEBHDAwvxTYPJPbVtXZVbWyqlYuWbJkpwuVJD24cYbHWmBFkgOTLKYb2mT1kNteBLw4yd5J9gZe3C+TJM2CsYVHVW0FTqb7pX898Omq2pDktCTHAiR5ZpIJ4NXAWUk29NveDvwpXQCtBU7bdvJckjR+zSfMd0VVrWHSifaqOnVgei1dl9RU254LnDvSAiVJQ/ExtJKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKajTU8khyd5IYkG5OsmmL9w5J8ql9/eZLl/fLlSe5KclX/9TfjrFuStL09xvVGSRYBZwIvAiaAtUlWV9V1A83eDNxRVU9IcgJwBvCaft13quqwcdUrSZreOI88jgA2VtWmqroHOB84blKb44BP9NMXAC9MkjHWKEkawjjDY3/gloH5iX7ZlG2qaivwQ2Dfft2BSdYn+WqSI0ddrCRpemPrtgKmOoKoIdvcCiyrqu8neQbw+SRPrqofbbdxchJwEsCyZctmoGRJ0lTGeeQxARwwML8U2DxdmyR7AI8Bbq+qu6vq+wBVdQXwHeDfT36Dqjq7qlZW1colS5aM4CNIkmC84bEWWJHkwCSLgROA1ZParAZO7KePB75cVZVkSX/CnSQHASuATWOqW5I0ydi6rapqa5KTgYuARcC5VbUhyWnAuqpaDXwc+GSSjcDtdAED8DzgtCRbgfuA36uq28dVuyRpe+M850FVrQHWTFp26sD0z4BXT7HdZ4HPjrxASdJQvMNcktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVKzsd4kKEkzafmqC++fvvH0Y2axkoXHIw9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1G2t4JDk6yQ1JNiZZNcX6hyX5VL/+8iTLB9a9p19+Q5KXjLNuSdL2xhYeSRYBZwIvBQ4BXpvkkEnN3gzcUVVPAP4SOKPf9hDgBODJwNHAx/rXkyTNgnEeeRwBbKyqTVV1D3A+cNykNscBn+inLwBemCT98vOr6u6q+jdgY/96kqRZMM7w2B+4ZWB+ol82ZZuq2gr8ENh3yG0lSWOyxxjfK1MsqyHbDLMtSU4CTupnf5LkhqYKH2g/4Hu7+BojlTNm/CXnzGcewWeb7jXnzGdutQs/o932M09niJ/F2D/zKP4NN9qZz/z4YRqNMzwmgAMG5pcCm6dpM5FkD+AxwO1DbktVnQ2cPVMFJ1lXVStn6vV2B37mhcHPvDCM8jOPs9tqLbAiyYFJFtOdAF89qc1q4MR++njgy1VV/fIT+quxDgRWAN8YU92SpEnGduRRVVuTnAxcBCwCzq2qDUlOA9ZV1Wrg48Ank2ykO+I4od92Q5JPA9cBW4G3V9V946pdkrS9cXZbUVVrgDWTlp06MP0z4NXTbPt+4P0jLfCBZqwLbDfiZ14Y/MwLw8g+c7peIUmShufwJJKkZobHFHY0jMp8kOSAJJcmuT7JhiTv6Jfvk+TiJN/uv+8927XOtCSLkqxP8k/9/IH9cDjf7ofHWTzbNc6kJHsluSDJt/r9/WvzfT8n+S/9v+trk5yX5OHzcT8nOTfJbUmuHVg25b5N5yP977VvJnn6rry34THJkMOozAdbgXdV1cHAs4G3959zFXBJVa0ALunn55t3ANcPzJ8B/GX/me+gGyZnPvkr4J+r6knA0+g++7zdz0n2B34fWFlVT6G7QOcE5ud+/lu6IZsGTbdvX0p3peoKuvvh/npX3tjweKBhhlHZ7VXVrVV1ZT/9Y7pfKPuz/RAxnwBeMTsVjkaSpcAxwDn9fIAX0A2HA/PsMyf5JeB5dFcyUlX3VNUPmOf7me5ioEf094s9EriVebifq+prdFemDppu3x4H/K/qXAbsleRxO/vehscDLbihUPrRiw8HLgf+XVXdCl3AAI+dvcpG4sPAHwA/7+f3BX7QD4cD829/HwRsAf5n31V3TpI9mcf7uar+L/BB4Ga60PghcAXzez8Pmm7fzujvNsPjgYYaCmW+SPIo4LPAO6vqR7NdzygleTlwW1VdMbh4iqbzaX/vATwd+OuqOhz4KfOoi2oqfR//ccCBwK8Ae9J12Uw2n/bzMGb037rh8UBDDYUyHyR5KF1w/H1Vfa5f/N1th7L999tmq74ReC5wbJIb6bojX0B3JLJX370B829/TwATVXV5P38BXZjM5/18FPBvVbWlqu4FPgc8h/m9nwdNt29n9Heb4fFAwwyjstvr+/o/DlxfVR8aWDU4RMyJwD+Ou7ZRqar3VNXSqlpOt1+/XFWvAy6lGw4H5t9n/n/ALUme2C96Id1IDfN2P9N1Vz07ySP7f+fbPvO83c+TTLdvVwO/3V919Wzgh9u6t3aGNwlOIcnL6P4i3TaMyrjvbB+5JL8O/AtwDb/o//9DuvMenwaW0f0nfHVVTT4ht9tL8nzg3VX18iQH0R2J7AOsB15fVXfPZn0zKclhdBcILAY2Ab9D94fjvN3PSf4EeA3dVYXrgbfQ9e/Pq/2c5Dzg+XSj534XeB/weabYt32QfpTu6qw7gd+pqnU7/d6GhySpld1WkqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4aMFKcl+Sq/pnPnwmySN34bX+OMm7++mvP0i7vZK8bWffZ4g6HpHkq/2jBUhyeJLTkyxO8rWB4TmkXWJ4aCG7q6oO65/5cA/we4Mr+2Ecmv+PVNVzHmT1XkBTeDTW8Sbgc1V1X1/L+qpa1T9e4BK6u66lXWZ4SJ1/AZ6QZHn/tL2PAVcCByR5fZJv9EcpZw38Vf/edE+c/BKwbewokvxkYPq3+6e2XZ3kk8DpwK/2r/WBvs0p/dHPtUne2S+bqo49k1zYv9a1SaYKgtcxMGZTkr/rh2KBbtiK183YT0wLmoewWvD6rpyXAv/cL3oi3bg/b0tyMN1f68+tqnv7X+avS7KBbnDFw+n+H11J98yIwdd9MvDeftvvJdkH+CXgKVV1WN/mGXRjTT2Lbsjsy5N8le5Jd/fX0bd9FbC5qo7p5x8z6f0WAwdV1Y0Di58KfLOfvhZ45s7/pKRf8MhDC9kjklwFrKMbQO7j/fKb+ietQTci6zOAtX3bF9I9YOlI4B+q6s7+OShTjbz8AuCCqvoewDQDD/56/zo/raqf0A0ffuQUdUA3iOVRSc5IcmRV/XDSa+0H/GDbTB8me257374r654kj97Bz0XaIY88tJDdte0IYJtu4FF+OrgI+ERVvWdSu3ey4wfpZMg20xmsg6r61/5I5WXAnyf5YlWdNtDkLuDhA/OHsP2z2gEeBvxsBzVJO+SRh/TgLgGOT/JYgCT7JHk88DXglf3VTY8Gfmuabf9jkn23bQv8GBj8y/9rwCv6Z0/sCbyS7vzLAyT5FeDOqvo7usesPn1wfVXdASxKsi1AngZcPbD9vsC2ByRJu8QjD+lBVNV1Sf478MX+iqd7gbdX1WVJPgVcBdzEFL/wq2pDkvcDX01yH7C+qt6Y5H8nuRb4QlX91yR/C3yj3+ycqlqf7rnykx0KfCDJz/s63jpFmy/SdYV9iS48Bru9fhNY0/gjkKbk8zykeSTJ4cApVfWGJFcCr6yqm/p1nwPeU1U3zGqRmhfstpLmkapaD1yWZD1w8UBwLAY+b3BopnjkIUlq5pGHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRm/x+GxQxOTbpwOAAAAABJRU5ErkJggg==\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(\"Predictors ($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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([13, 31, 38, 55, 77, 89]),)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "actual_model = data.mean(axis=0) > 0\n",
    "actual_model.nonzero()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the only predictors with non zero posterior are the ones belonging to our original boolean function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "models = pd.DataFrame(data).assign(count=-1).groupby(list(range(100))).count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAENCAYAAAD34uk0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFuRJREFUeJzt3X+0ZWV93/H3J4P4uwZlrIYBByIqqATiiGksGhUUgwGT6nKoukglndaKmhKbjLVBO1lJ8UdN0ooVFk5xGQsi+GNWHYuIRLuWIjMCAgNSR0LgCg2jQ9QoigPf/rH30DPXe4fzMHP2Ocy8X2vddc9+9rPv+c6Pez93P3vv50lVIUlSi1+YdgGSpIcew0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUrN9pl3ApOy///61fPnyaZchSQ8pX//6179bVUsfqN8eGx7Lly9n48aN0y5Dkh5SkvztOP0ctpIkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ122OfMNfwlq/+7IM+9pYzT9iNlUiaNM88JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1GzQ8khyf5KYkm5OsXmD/6UluSHJtksuSPGVk371Jruk/1g1ZtyRpR4M9YZ5kCXAWcBwwB2xIsq6qbhjpdjWwoqp+nOSNwHuA1/T77q6qI4eqV5K0uCHPPI4GNlfVzVV1D3ABcNJoh6q6vKp+3G9eASwbsD5J0piGDI8DgNtGtuf6tsWcCnxuZPsRSTYmuSLJKxc6IMmqvs/GLVu27HrFkqQFDTkxYhZoqwU7Jq8DVgAvHGk+qKpuT3II8MUk11XVt3f4YlXnAOcArFixYsGvLUnadUOeecwBB45sLwNun98pybHAO4ATq+qn29ur6vb+883AXwNHTbJYSdLihgyPDcChSQ5Osi+wEtjhrqkkRwFn0wXHnSPt+yV5eP96f+D5wOiFdknSgAYbtqqqbUlOAy4BlgBrq2pTkjXAxqpaB7wXeAzwiSQAt1bVicBhwNlJ7qMLvDPn3aUlSRrQoItBVdV6YP28tjNGXh+7yHFfAZ492eokSePyCXNJUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1GzQ8EhyfJKbkmxOsnqB/acnuSHJtUkuS/KUkX2nJPlW/3HKkHVLknY0WHgkWQKcBbwcOBw4Ocnh87pdDayoqiOAi4D39Mc+Hngn8DzgaOCdSfYbqnZJ0o6GPPM4GthcVTdX1T3ABcBJox2q6vKq+nG/eQWwrH/9MuDSqtpaVXcBlwLHD1S3JGmeIcPjAOC2ke25vm0xpwKfazk2yaokG5Ns3LJlyy6WK0lazJDhkQXaasGOyeuAFcB7W46tqnOqakVVrVi6dOmDLlSStHNDhscccODI9jLg9vmdkhwLvAM4sap+2nKsJGkYQ4bHBuDQJAcn2RdYCawb7ZDkKOBsuuC4c2TXJcBLk+zXXyh/ad8mSZqCfYZ6o6raluQ0uh/6S4C1VbUpyRpgY1WtoxumegzwiSQAt1bViVW1Ncmf0AUQwJqq2jpU7ZKkHQ0WHgBVtR5YP6/tjJHXx+7k2LXA2slVJ0kal0+YS5KaGR6SpGaDDltJemhYvvqzu3T8LWeesJsq0azyzEOS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1aw6PJI9OsmQSxUiSHhoeMDyS/EKSf57ks0nuBL4J3JFkU5L3Jjl08mVKkmbJOGcelwO/DLwdeFJVHVhVTwSOAa4AzkzyugnWKEmaMeOsJHhsVf1sfmNVbQUuBi5O8rDdXpkkaWaNc+bxS4vtSLICYKFwkSTtucYJj01J/mj0InmSxyb5APDVyZUmSZpV44THHPCfgGuSHJPk1XQXzf8N8JVJFidJmk3jhMezgXfRXTT/a+ACYF/g96rqhROrTJI0sx4wPPrrGZ8BNgPpP64GLp1saZKkWfWAd1sleR/wFqCA9wFPBl4L3JDkXVX1nydb4nQsX/3ZB33sLWeesBsrkaTZM86w1enAlcBzquoPq+r1wHHAHcB7JlmcJGk2jfOcx6qqOne0oaouS/Is4N9PpixJ0iwb55rHuYu031NV72p5syTHJ7kpyeYkqxfY/4IkVyXZluRV8/bdm+Sa/mNdy/tKknavcc48dov+OZGz6Ia85oANSdZV1Q0j3W4Ffhd42wJf4u6qOnLihUqSHtCDnpI9yZOTPLzhkKOBzVV1c1XdQ3fL70mjHarqlqq6FrjvwdYlSZq8XVnP46PAN/u7scZxAHDbyPZc3zauRyTZmOSKJK9sOE6StJuNc6vuI6vq7vntVXVskgCHj/leWaCtxjwW4KCquj3JIcAXk1xXVd+eV+sqYBXAQQcd1PClJUktxjnz+GaS31loR3U2jflec8CBI9vLgNvHPJaqur3/fDPdk+5HLdDnnKpaUVUrli5dOu6XliQ1Gic8ngR8Isnnkjx1F95rA3BokoOT7AusBMa6ayrJftuvryTZH3g+cMPOj5IkTco4d1s9E3g/8ArguiQfYuTaRVW9f5w3qqptSU4DLgGWAGuralOSNcDGqlqX5LnAp4D9gN9K8h+r6pnAYcDZSe6jC7wz592lJUka0AOGR1VtBk5MchbwRuCto7vpgmUsVbUeWD+v7YyR1xvohrPmH/cVugkaJUkzYJwL5kcCf0G37OwP6ea3unXCdUmSZtg4w1Yb6YaKPg6cXlV3TLYkSdKsGyc8bgJOq6rLJ12MJOmhYZzw+JWq2jbxSiRJDxnjTIxocEiSdrAr05NIkvZShockqdk4t+qeDzxlZ13oZir59d1WlRa1K8vjgkvkSto9xnlI8OQhCpEkPXQ4bCVJamZ4SJKaGR6SpGYPGB5JXpzExTEkSfcb5wnzLwB39tOhXw9cC1zXf76hqn46wfokSTNonPB4C/AG4ELgK8DTgecAv0u3zsaTJlWcJGk2jTM9yQfoVu4ruqnZfwa8tapeVFUGhyTthca6YF5Vd1fVu4HfAJ4KXJnkeZMsTJI0u8Z5wvwYuuGpZ/Sfn0i3KNQTJluaJGlWjXPN40vAN4Dzgf9SVbdMtCJJ0swbJzzeSLd++AnAHyT5Ht3dVtcB11fVpydYnyRpBo0THudUVW3fSLIMOIIuUP4Z8OkkGe0jSdqzjRMelye5GPhMVd1aVXPAXJIvAMck+QhwOXDeBOuUJM2QccLjeLrnPM5PcghwF/AIYAnweeDPq+qayZUoSZo140zJ/hPgg8AHkzwM2B+4u6r+ftLFSZJm0zhzW52S5LtJtgLnAv9gcEjS3m2cYas/Bo4DvgO8Gfiz/rOkXbQrK0O6KqSmaZzw+EFVXd2//uMkX5tkQZKk2TdOeDw5ySrgRuCbwMMmW5IkadaNEx7vpHuu47V0z3Y8Jsl6uqfOr62q8ydYnyRpBo1zt9U5o9vzHhL8TbppSyRJe5HmZWiraq6q1lfVu6vq9S3HJjk+yU1JNidZvcD+FyS5Ksm2JK+at++UJN/qP05prVuStPsMtoZ5kiXAWcDLgcOBk5McPq/brXSLTP2Pecc+nm747HnA0cA7k+w36ZolSQsbLDzofuhvrqqbq+oe4ALgpNEOVXVLVV0L3Dfv2JcBl1bV1qq6C7iU7sl3SdIUDBkeBwC3jWzP9W2TPlaStJsNGR5ZoG3cmXjHOjbJqiQbk2zcsmVLU3GSpPENGR5zwIEj28uA23fnsVV1TlWtqKoVS5cufdCFSpJ2bsjw2AAcmuTgJPsCK4F1Yx57CfDSJPv1F8pf2rdJkqZgsPCoqm3AaXQ/9G8ELqyqTUnWJDkRIMlzk8wBrwbOTrKpP3Yr8Cd0AbQBWNO3SZKmYJwnzHebqloPrJ/XdsbI6w10Q1ILHbsWWDvRAiVJYxly2EqStIcwPCRJzQwPSVIzw0OS1MzwkCQ1G/RuK2kadmWpV3C5V2khnnlIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqNmh4JDk+yU1JNidZvcD+hyf5eL//a0mW9+3Lk9yd5Jr+40ND1i1J2tE+Q71RkiXAWcBxwBywIcm6qrphpNupwF1V9dQkK4F3A6/p9327qo4cql5J0uKGPPM4GthcVTdX1T3ABcBJ8/qcBHykf30R8JIkGbBGSdIYhgyPA4DbRrbn+rYF+1TVNuD7wBP6fQcnuTrJl5IcM+liJUmLG2zYCljoDKLG7HMHcFBVfS/Jc4BPJ3lmVf1gh4OTVcAqgIMOOmg3lCxJWsiQZx5zwIEj28uA2xfrk2Qf4HHA1qr6aVV9D6Cqvg58G3ja/DeoqnOqakVVrVi6dOkE/giSJBg2PDYAhyY5OMm+wEpg3bw+64BT+tevAr5YVZVkaX/BnSSHAIcCNw9UtyRpnsGGrapqW5LTgEuAJcDaqtqUZA2wsarWAR8GPppkM7CVLmAAXgCsSbINuBf411W1dajaJUk7GvKaB1W1Hlg/r+2Mkdc/AV69wHEXAxdPvEBJ0lh8wlyS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNBl3DfG+1fPVnH/Sxt5x5wm6sRJJ2D888JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ18yFBSdoNduVhYHjoPRBseOzF9rb/7Hu6vWUmg73lzznrBh22SnJ8kpuSbE6yeoH9D0/y8X7/15IsH9n39r79piQvG7JuSdKOBguPJEuAs4CXA4cDJyc5fF63U4G7quqpwJ8D7+6PPRxYCTwTOB74YP/1JElTMOSw1dHA5qq6GSDJBcBJwA0jfU4C3tW/vgj4QJL07RdU1U+Bv0myuf96Xx2odkm7YFaHmhy6ffCGDI8DgNtGtueA5y3Wp6q2Jfk+8IS+/Yp5xx4wuVI1bbP6w0ZSJ1U1zBslrwZeVlW/12+/Hji6qt480mdT32eu3/423RnGGuCrVfVXffuHgfVVdfG891gFrOo3nw7cNLJ7f+C7k/iz7QazWtus1gWzW9us1gWzW9us1gWzW9sk63pKVS19oE5DnnnMAQeObC8Dbl+kz1ySfYDHAVvHPJaqOgc4Z6E3T7KxqlY86OonaFZrm9W6YHZrm9W6YHZrm9W6YHZrm4W6hrzbagNwaJKDk+xLdwF83bw+64BT+tevAr5Y3anROmBlfzfWwcChwJUD1S1JmmewM4/+GsZpwCXAEmBtVW1KsgbYWFXrgA8DH+0viG+lCxj6fhfSXVzfBrypqu4dqnZJ0o4GfUiwqtYD6+e1nTHy+ifAqxc59k+BP92Ft19wOGtGzGpts1oXzG5ts1oXzG5ts1oXzG5tU69rsAvmkqQ9hxMjSpKa7RXh8UDTokxDkgOTXJ7kxiSbkrx12jWNSrIkydVJ/ue0axmV5BeTXJTkm/3f3T+Zdk3bJfm3/b/l9UnOT/KIKdayNsmdSa4faXt8kkuTfKv/vN+M1PXe/t/z2iSfSvKLQ9e1WG0j+96WpJLsPyt1JXlz/3NtU5L3DF3XHh8eY06LMg3bgD+oqsOAXwPeNCN1bfdW4MZpF7GAvwT+V1U9A/gVZqTGJAcAbwFWVNWz6G4KWTnFks6jm8pn1Grgsqo6FLis3x7aefx8XZcCz6qqI4D/A7x96KJ65/HztZHkQOA44NahC+qdx7y6kryIbuaNI6rqmcD7hi5qjw8PRqZFqap7gO3TokxVVd1RVVf1r39I90NwJp6aT7IMOAE4d9q1jEryj4AX0N2VR1XdU1V/P92qdrAP8Mj+GaVHscCzSEOpqi/T3bE46iTgI/3rjwCvHLQoFq6rqj5fVdv6zSvonuMa3CJ/Z9DNs/eHwFQuEC9S1xuBM/spm6iqO4eua28Ij4WmRZmJH9Lb9bMHHwV8bbqV3O8v6L5Z7pt2IfMcAmwB/ns/pHZukkdPuyiAqvoO3W9/twJ3AN+vqs9Pt6qf84+r6g7ofnkBnjjlehbyBuBz0y5iuyQnAt+pqm9Mu5Z5ngYc088+/qUkzx26gL0hPLJA28zcYpbkMcDFwO9X1Q9moJ5XAHdW1denXcsC9gF+FfhvVXUU8COmM/Tyc/rrBycBBwO/BDw6yeumW9VDS5J30A3nfmzatQAkeRTwDuCMB+o7BfsA+9ENef874MJ+EtnB7A3hMdbUJtOQ5GF0wfGxqvrktOvpPR84McktdEN8L07yV9Mt6X5zwFxVbT9Du4guTGbBscDfVNWWqvoZ8Eng16dc03x/l+TJAP3nwYc6FpPkFOAVwGtrdp4f+GW6Xwa+0X8/LAOuSvKkqVbVmQM+WZ0r6UYJBr2YvzeExzjTogyu/y3hw8CNVfX+adezXVW9vaqWVdVyur+rL1bVTPwGXVX/F7gtydP7ppew45T+03Qr8GtJHtX/276EGbmYP2J0+p9TgM9MsZb7JTke+CPgxKr68bTr2a6qrquqJ1bV8v77YQ741f7/4bR9GngxQJKnAfsy8ASOe3x49Bfitk+LciNwYVVtmm5VQPcb/uvpfrO/pv/4zWkX9RDwZuBjSa4FjgT+bMr1ANCfDV0EXAVcR/e9NbWngJOcT7fezdOTzCU5FTgTOC7Jt+juHjpzRur6APBY4NL+++BDQ9e1k9qmbpG61gKH9LfvXgCcMvQZm0+YS5Ka7fFnHpKk3c/wkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9pAf3aDR8d2d4nyZaW9U2SvCvJ28bo9w+LtD+yn/RuSb/9r/q6XjjS57S+7aVJvtzP6CtNnOEhLexHwLOSPLLfPg74zsA1vIFu/qJ7++0jgGuBw+D+iftOpZtp+Bt0a3S8ZuAatZcyPKTFfY5uXROAk4Hzt+9Icnq/YuD1SX5/pP0d/epuXwCePvrFkrwuyZX9FBxnbz+j2InXsuP8U8/ua3hGv/0W4BPAfVX1d3TzHb22/Y8ptTM8pMVdAKxMt5zsEfTrrSR5DvAvgOfRTYn9L5Mc1bevpFub5XeA+9dYSHIY3VnB86vqSOBedvKDvp/E85CqumWk+TDgQuAZSR7Xf72vANuXJ71+9D2lSXJ8VFpEVV3bL9R1MrB+ZNc/BT5VVT8CSPJJ4Bi6X8Y+tX1m2CSjsze/BHgOsKFfduGR7HxK9P2B+1dJ7JdC/V5V3ZzkiXSLdf1XukWBru3rvTfJPUke269OKU2M4SHt3Dq6FQJ/A3hC37azRXcWm2k0wEeqatz1ue8GHjGyfQTdbL0AP6Rb0/poulUfrxrp93DgJ2O+h/SgOWwl7dxaYE1VXTfS9mXglf3aHY8Gfhv43337b/d3ST0W+K2RYy4DXtWfNZDk8UmestibVtVdwJJ+yAy66x3ba3gvcFp/If3Z9GceSZ4AbF+MSpoozzyknaiqOeAv57VdleQ84Mq+6dyquhogyceBa4C/pQuU7cfckOQ/AJ9P8gvAz4A39f0W83m6IbIv0IXExf3XGr1d+HD+/4JYL2LH4TVpYlzPQ5pRSY4CTq+q14/Z/5PA26vqpslWJjlsJc2s/mzm8jFu6d1+d9anDQ4NxTMPSVIzzzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LU7P8B2MW/miuGResAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "model_counts = models.loc[:, 'count'].values\n",
    "ax.bar(np.arange(1, model_counts.size + 1), model_counts / model_counts.sum())\n",
    "ax.set_xlabel(\"Model ($M$)\")\n",
    "ax.set_ylabel(r\"$P(M \\mid \\mathbf{Y}, \\mathbf{X})$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like the 5th model has a very high posterior. Let's see which predictors it contains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([38]),)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_estimate = models.reset_index().iloc[4,:100].values > 0\n",
    "model_estimate.nonzero()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is disappointing. Let's compare the log posterior of actual and estimated models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-647.3964116216055\n",
      "-734.0178634863765\n"
     ]
    }
   ],
   "source": [
    "print(main.posterior(Y, X, model_estimate, 0.1, 4, np.e, 4))\n",
    "print(main.posterior(Y, X, actual_model, 0.1, 4, np.e, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The actual model does have a lower posterior."
   ]
  }
 ],
 "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
}