{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2 Markov Chain Monte Carlo\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://docs.pymc.io/notebooks/GLM-logistic.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import pymc3 as pm\n",
    "from sklearn import preprocessing\n",
    "from scipy import stats\n",
    "from scipy import special\n",
    "\n",
    "# Use a theano shared variable to be able to exchange the data the model runs on\n",
    "from theano import shared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Recap\n",
    "---\n",
    "\n",
    "\\begin{equation}\n",
    "P(\\boldsymbol{\\theta}\\ |\\ D) =  \\frac{P(D\\ |\\ \\boldsymbol{\\theta})\\cdot P(\\boldsymbol{\\theta})}{P(D)}\n",
    "\\end{equation}\n",
    "\n",
    "---\n",
    "### Example\n",
    "As an introductory example, we observe two weights, 65kg, and 75kg.\n",
    "\n",
    "We would like to learn something about the distribution of the weights.\n",
    "\n",
    "![](images/weights.png)\n",
    "\n",
    "---\n",
    "How do we model this?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 The data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = np.array([65, 75])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 Parameter space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets creat a grid of our two parameters\n",
    "mu = np.linspace(60, 90)\n",
    "sigma = np.linspace(0, 10)[::-1]\n",
    "\n",
    "mm, ss = np.meshgrid(mu, sigma)  # just broadcasted parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 Likelihood\n",
    "$$ \\prod_{i=0}^N{P(D_i|\\theta)} $$\n",
    "$$ P(D_i|\\theta) \\sim N(\\theta) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\sigma$')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 9652 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "likelihood = stats.norm(mm, ss).pdf(weights[0]) * stats.norm(mm, ss).pdf(weights[1])\n",
    "aspect = mm.max() / ss.max() / 3\n",
    "extent = [mm.min(), mm.max(), ss.min(), ss.max()]\n",
    "# extent = left right bottom top\n",
    "plt.imshow(likelihood, cmap='Reds', aspect=aspect, extent=extent) \n",
    "plt.xlabel(r'$\\mu$')\n",
    "plt.ylabel(r'$\\sigma$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4 Prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\sigma$')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 7868 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prior = stats.norm(80, 3).pdf(mm) * stats.cauchy(0, 2).pdf(ss)\n",
    "\n",
    "plt.imshow(prior, cmap='Greens', aspect=aspect, extent=extent)\n",
    "plt.xlabel(r'$\\mu$')\n",
    "plt.ylabel(r'$\\sigma$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.5 Posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\sigma$')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 8160 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "unnormalized_posterior = prior * likelihood\n",
    "posterior = unnormalized_posterior / np.nan_to_num(unnormalized_posterior).sum()\n",
    "plt.imshow(posterior, cmap='Blues', aspect=aspect, extent=extent)\n",
    "plt.xlabel(r'$\\mu$')\n",
    "plt.ylabel(r'$\\sigma$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.6 We can do the same with pymc3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([65, 75])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sequential sampling (1 chains in 1 job)\n",
      "CompoundStep\n",
      ">Metropolis: [sigma]\n",
      ">Metropolis: [mu]\n",
      "100%|██████████| 10500/10500 [00:02<00:00, 4563.95it/s]\n",
      "Only one chain was sampled, this makes it impossible to run some convergence checks\n"
     ]
    }
   ],
   "source": [
    "with pm.Model():\n",
    "    mu = pm.Normal('mu', mu=80, sd=3)\n",
    "    sigma = pm.HalfCauchy('sigma', 2)\n",
    "\n",
    "    observed = pm.Normal('observed', mu=mu, sd=sigma, observed=weights)\n",
    "    trace = pm.sample(draws=10000, chains=1, step=pm.Metropolis())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 47040 characters]\n",
      "text/plain": [
       "<Figure size 864x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot(trace);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 10.0)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 42712 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(2, sharex=True, sharey=True, figsize = (8, 8))\n",
    "axes[0].imshow(posterior, cmap='Blues', extent=extent, aspect=4)\n",
    "axes[1].scatter(trace['mu'], trace['sigma'], alpha=0.01)\n",
    "plt.xlim(0, mm.max())\n",
    "plt.ylim(0, ss.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.7 Conclusions\n",
    "\n",
    "- Posterior is proportional to likelihood multiplied with prior.\n",
    "- We 'showed' (with plots) that samples obtained with `pymc3` are from this distribution.\n",
    "- `pymc3` can be used to approximate any posterior, even if we can not validate this with plots or math."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Intermezzo: Logistic regression \n",
    "\n",
    "Let's do a Baysian Logistic regression. What is the chance I'll make \\$50k?\n",
    "\n",
    "First we load and preprocess some data. We use **age, education and working hours** as features $X$ and we try to predict:\n",
    "\n",
    "$$ y = \\begin{cases}\n",
    "      1, & \\text{if income > \\$50K}  \\\\\n",
    "      0, & \\text{otherwise}\n",
    "    \\end{cases}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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>age</th>\n",
       "      <th>educ</th>\n",
       "      <th>hours</th>\n",
       "      <th>more_than_50k</th>\n",
       "      <th>age_squared</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>39</td>\n",
       "      <td>13</td>\n",
       "      <td>40</td>\n",
       "      <td>0</td>\n",
       "      <td>1521</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>50</td>\n",
       "      <td>13</td>\n",
       "      <td>13</td>\n",
       "      <td>0</td>\n",
       "      <td>2500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>38</td>\n",
       "      <td>9</td>\n",
       "      <td>40</td>\n",
       "      <td>0</td>\n",
       "      <td>1444</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>53</td>\n",
       "      <td>7</td>\n",
       "      <td>40</td>\n",
       "      <td>0</td>\n",
       "      <td>2809</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>37</td>\n",
       "      <td>14</td>\n",
       "      <td>40</td>\n",
       "      <td>0</td>\n",
       "      <td>1369</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   age  educ  hours  more_than_50k  age_squared\n",
       "0   39    13     40              0         1521\n",
       "1   50    13     13              0         2500\n",
       "2   38     9     40              0         1444\n",
       "3   53     7     40              0         2809\n",
       "5   37    14     40              0         1369"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv(\"https://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data\", header=None, names=['age', 'workclass', 'fnlwgt',\n",
    "                'education-categorical', 'educ',\n",
    "                'marital-status', 'occupation',\n",
    "                'relationship', 'race', 'sex',\n",
    "                'captial-gain', 'capital-loss',\n",
    "                'hours', 'native-country',\n",
    "                'income'])\n",
    "data = data[data['native-country'] == ' United-States']\n",
    "data = data[~pd.isnull(data['income'])]\n",
    "data = data.assign(more_than_50k=(data['income'] == ' >50K') * 1)\n",
    "education_map = dict(zip(list(map(lambda x: x.strip(), data['education-categorical'].unique())), data['educ'].unique()))\n",
    "data = data[['age', 'educ', 'hours', 'more_than_50k']]\n",
    "data = data.assign(age_squared=data['age']**2)\n",
    "\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a correlation table for the features and the labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fe402eed1d0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 14740 characters]\n",
      "text/plain": [
       "<Figure size 1152x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(16, 8))\n",
    "corr = data.corr()\n",
    "mask = np.zeros_like(corr, dtype=np.bool)\n",
    "# indices for upper triangle of matrix\n",
    "idx = np.triu_indices_from(corr)\n",
    "mask[idx] = True\n",
    "sns.heatmap(corr, mask=mask, vmax=.3, cmap=sns.color_palette('Reds'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Model\n",
    "Now we need to think of a model that can explain this data. It's often easiest to start with the labels. The labels are $y \\in \\{0, 1\\}$. \n",
    "\n",
    "## Exercise\n",
    "* What is a distribution that can generate these labels?\n",
    "* What is/ are the parameter(s) of this distribution?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise answers\n",
    "\n",
    "A Bernoulli distribution represents binary outxomes. The Bernoulli has a parameter $p$, where \n",
    "\n",
    "$p = P(x=1)$\n",
    "\n",
    "* What is the formulation of our model?\n",
    "\n",
    "Hint: we need a link function $\\frac{1}{1 + e^{-z}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ y \\sim Bernoulli(p) $$\n",
    "$$ p = \\frac{1}{1 + e^{-z}} $$\n",
    "$$ z = \\beta X $$\n",
    "$$ \\beta \\sim N(0, 15) $$\n",
    "\n",
    "* Now we need to set $P(\\theta)$ e.g. priors on $\\beta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe407bd1a20>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 15228 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-75, 75)\n",
    "plt.plot(x, stats.norm(0, 15).pdf(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:595: DataConversionWarning: Data with input dtype int64 was converted to float64 by MinMaxScaler.\n",
      "  warnings.warn(msg, DataConversionWarning)\n"
     ]
    }
   ],
   "source": [
    "n = 2000\n",
    "x = data[['age', 'educ', 'hours', 'age_squared']].values[:n]\n",
    "y = data['more_than_50k'].values[:n]\n",
    "\n",
    "scaler = preprocessing.MinMaxScaler()\n",
    "scaled_x = scaler.fit_transform(x)\n",
    "\n",
    "# add one for intercept\n",
    "scaled_x = np.hstack((np.ones(x.shape[0])[:, None], scaled_x))\n",
    "shared_x = shared(scaled_x) # This make it possible to interact date with the Theano computation graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shared_x.get_value().max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    # 3 features and an intercept is 4 dimensions\n",
    "    beta = pm.Normal('beta', mu=0, sd=15, shape=5)\n",
    "    observed = pm.Bernoulli('observed', p=pm.math.invlogit(pm.math.dot(beta, shared_x.T)), observed=y)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Have we got an uninformative prior?\n",
    "Below we sample predictions from our posterior. We are uninformative if $P(>50K) \\sim 0.5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    prior_pc = pm.sample_prior_predictive()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fe4053925c0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 13808 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(prior_pc['observed'].mean(0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta]\n",
      "Sampling 4 chains: 100%|██████████| 4000/4000 [00:19<00:00, 200.58draws/s]\n",
      "The acceptance probability does not match the target. It is 0.8927220632971556, but should be close to 0.8. Try to increase the number of tuning steps.\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    trace = pm.sample()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7fe3f47d2320>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fe3f47f2278>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 64408 characters]\n",
      "text/plain": [
       "<Figure size 864x144 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([<matplotlib.axes._subplots.AxesSubplot object at 0x7fe3f4707780>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x7fe3f46cac50>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x7fe3f46f3e48>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x7fe3f46a00b8>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x7fe3f46492e8>],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 36568 characters]\n",
      "text/plain": [
       "<Figure size 864x540 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#  intercept, 'age', 'educ', 'hours', 'age_squared'\n",
    "pm.plot_posterior(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe3f44d52e8>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 17916 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "trace_age = trace['beta'][:, 1]\n",
    "trace_education = trace['beta'][:, 2]\n",
    "trace_hours = trace['beta'][:, 3]\n",
    "trace_age_squared = trace['beta'][:, 4]\n",
    "sns.distplot(trace_age_squared, label='age^2')\n",
    "sns.distplot(trace_age, label='age')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.4 Examine the the joint distribution of Betas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(336.9714285714286, 0.5, 'beta_hours')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 20704 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.jointplot(trace_age, trace_hours, kind='kde')\n",
    "plt.xlabel(\"beta_age\")\n",
    "plt.ylabel(\"beta_hours\")\n",
    "\n",
    "# and so forth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Posterior Predictive Check\n",
    "With the posterior predictive check we sample predicitons for our inferred model with the posterior distribution for $\\theta$.\n",
    "Below we make a plot of our probability of making more than 50K. $P(>50K| m)$. Note that this is the probability given our inferred model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:2: DeprecationWarning: sample_ppc() is deprecated.  Please use sample_posterior_predictive()\n",
      "  \n",
      "100%|██████████| 1000/1000 [00:00<00:00, 1446.93it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1000, 2000)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with model:\n",
    "    ppc = pm.sample_ppc(trace, samples=1000)['observed']\n",
    "ppc.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.251, 3, 'Average P(>50K)')"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 18300 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title('Overall P(>50K)')\n",
    "sns.distplot(ppc.mean(1))\n",
    "plt.vlines(y.sum() / y.shape[0], 0, 35)\n",
    "plt.annotate('Average P(>50K)', (y.sum() / y.shape[0] + 0.002, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.04109589, 0.8       , 0.39795918, 0.01421073],\n",
       "       [1.        , 0.04109589, 0.8       , 0.12244898, 0.01421073],\n",
       "       [1.        , 0.04109589, 0.53333333, 0.39795918, 0.01421073],\n",
       "       ...,\n",
       "       [1.        , 0.04109589, 0.13333333, 0.29591837, 0.01421073],\n",
       "       [1.        , 0.04109589, 0.53333333, 0.39795918, 0.01421073],\n",
       "       [1.        , 0.04109589, 0.6       , 0.47959184, 0.01421073]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def prepare_shared_x(col_idx, new_value):\n",
    "    \"\"\"\n",
    "    :param col_idx: (int) 0 = age, 1 = education, 2 = hours\n",
    "    :param new_value: (flt) Replace the column values with this new value.\n",
    "    \"\"\"\n",
    "    shared_x.set_value(scaled_x)\n",
    "    feat = np.array(x)\n",
    "    feat[:, col_idx] = new_value\n",
    "    if col_idx == 0: # age is also squared\n",
    "        feat[:, -1] = new_value**2\n",
    "#     feat[:, 1] = 5\n",
    "    feat = scaler.transform(feat)\n",
    "    feat = np.hstack((np.ones(feat.shape[0])[:, None], feat))\n",
    "    shared_x.set_value(feat)\n",
    "    \n",
    "prepare_shared_x(0, 20)\n",
    "shared_x.get_value()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "quants = []\n",
    "age = np.arange(15, 80, step=5)\n",
    "with model:\n",
    "    for a in age:\n",
    "        prepare_shared_x(0, a)\n",
    "        ppc = pm.sample_posterior_predictive(trace, samples=100, progressbar=False)['observed']\n",
    "        quants.append(np.quantile(ppc.mean(1), [0.05, 0.5, 0.95]))\n",
    "quants = np.stack(quants)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x7fe3f43c86d8>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 21052 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(age, quants[:, 1])\n",
    "plt.ylabel('P(>50k)')\n",
    "plt.xlabel('age')\n",
    "plt.fill_between(age, quants[:, 0], quants[:, 2], alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Bachelors': 13,\n",
       " 'HS-grad': 9,\n",
       " '11th': 7,\n",
       " 'Masters': 14,\n",
       " 'Some-college': 10,\n",
       " 'Assoc-acdm': 12,\n",
       " 'Doctorate': 16,\n",
       " '9th': 5,\n",
       " 'Assoc-voc': 11,\n",
       " '10th': 6,\n",
       " '7th-8th': 4,\n",
       " 'Prof-school': 15,\n",
       " '1st-4th': 2,\n",
       " 'Preschool': 1,\n",
       " '5th-6th': 3,\n",
       " '12th': 8}"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "education_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "* Make a probability plot hours worked and P(>50K).\n",
    "* Plot P(>50K) for HS-grad, Bachelors and Masters.\n",
    "\n",
    "Note, columns are: 0 = age, 1 = education, 2 = hours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x7fe3f39801d0>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 16584 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "quants = []\n",
    "hours = np.arange(1, 80, step=5)\n",
    "with model:\n",
    "    for a in hours:\n",
    "        prepare_shared_x(2, a)\n",
    "        ppc = pm.sample_posterior_predictive(trace, samples=100, progressbar=False)['observed']\n",
    "        quants.append(np.quantile(ppc.mean(1), [0.05, 0.5, 0.95]))\n",
    "quants = np.stack(quants)\n",
    "plt.plot(age, quants[:, 1])\n",
    "plt.ylabel('P(>50k)')\n",
    "plt.xlabel('hours')\n",
    "plt.fill_between(age, quants[:, 0], quants[:, 2], alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe3f3c77780>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 43720 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "age = np.arange(1, 80, step=5)\n",
    "with model:\n",
    "    for education in ['HS-grad', 'Bachelors', 'Masters']:\n",
    "        quants = []\n",
    "        for a in age:\n",
    "            # prepare shared\n",
    "            shared_x.set_value(scaled_x)\n",
    "            feat = np.array(x)\n",
    "            feat[:, 0] = a\n",
    "            feat[:, -1] = a**2\n",
    "            feat[:, 1] = education_map[education]\n",
    "            feat = scaler.transform(feat)\n",
    "            feat = np.hstack((np.ones(feat.shape[0])[:, None], feat))\n",
    "            shared_x.set_value(feat)\n",
    "            ppc = pm.sample_posterior_predictive(trace, samples=100, progressbar=False)['observed']\n",
    "            quants.append(np.quantile(ppc.mean(1), [0.05, 0.5, 0.95]))\n",
    "            \n",
    "        quants = np.stack(quants)\n",
    "        plt.plot(age, quants[:, 1], label=education)\n",
    "        plt.fill_between(age, quants[:, 0], quants[:, 2], alpha=0.5)\n",
    "        \n",
    "plt.ylabel('P(>50k)')\n",
    "plt.xlabel('age')\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Markov Chain Monte Carlo [MCMC]\n",
    "\n",
    "* Markov Chain: A stochastic model describing a sequence of possible events in which the probability of each event depends only on the state attained in the previous event.\n",
    "* Monte Carlo: Probablistic simulation of events with stochastic variables.\n",
    "\n",
    "\n",
    "Or as I like to put it:\n",
    "\n",
    "*\"\n",
    "Markov chain Monte Carlo (MCMC) is a technique for estimating by simulation the expectation of a statistic in a complex model. Successive random selections form a Markov chain, the stationary distribution of which is the target distribution. It is particularly useful for the evaluation of posterior distributions in complex Bayesian models. In the Metropolis–Hastings algorithm, items are selected from an arbitrary “proposal” distribution and are retained or not according to an acceptance rule. The Gibbs sampler is a special case in which the proposal distributions are conditional distributions of single components of a vector parameter. Various special cases and applications are considered.\n",
    "\"* (Encyclopedia of Biostatistics)\n",
    "\n",
    "![](images/donotunderstand.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 King Markov\n",
    "[statistical rethinking](http://web.stanford.edu/class/stats305a/MCMC.html)\n",
    "Let's take a look at an analogy to understand what is really going on.\n",
    "\n",
    "![](images/islands.png)\n",
    "\n",
    "* King Markov has many islands and amnesia. He doesn't know how many islands he has nor the number of residents living on those islands.\n",
    "* King Markov wants to visit islands proportional to number of people living on it. For example, with two islands, with distribution `[1, 2]`, spend twice as much time on the second island as one the first.\n",
    "* He doesn't recognize islands from earlier visits.\n",
    "* He can only count the current island and one other island.\n",
    "\n",
    "**GOAL:** \n",
    "Find the distribution of number of residents on each island.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1.1 Simple algorithm\n",
    "* Count the number of inhabitants, stay that amount of day, then move to the next one.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 7676 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fe3f42aa7b8>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 8592 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "distribution = np.array([100, 200, 30, 10, 200, 50, 60, 80, 10, 1])\n",
    "\n",
    "def generate_path(distribution, start=0):\n",
    "    current_position = start\n",
    "    while True:\n",
    "        for _ in range(distribution[current_position]):  # stay the number of days\n",
    "            yield current_position  # return the current island to the distribution stack\n",
    "        current_position = (current_position + 1) % len(distribution)\n",
    "\n",
    "gen = generate_path(distribution, 0)\n",
    "path = np.array([next(gen) for _ in range(1200)])\n",
    "\n",
    "plt.plot(path)\n",
    "plt.show()\n",
    "(pd.Series(path).value_counts(normalize=1)\n",
    " .reindex(np.arange(len(distribution)))\n",
    " .rename('approx').to_frame()\n",
    " .assign(true=distribution / distribution.sum())\n",
    ").plot.bar()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After 1200 iterations we have visited al the islands, but the approximation is not very good.\n",
    "Below we show a plot of a situation with three islands and the approximations we get after **n** days of traveling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe3f41e1668>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 56088 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "distribution = np.array([50, 100, 10])\n",
    "gen = generate_path(distribution, 0)\n",
    "\n",
    "n = 1000\n",
    "path = np.array([next(gen) for _ in range(n)])\n",
    "\n",
    "plt.figure(figsize=(16, 6))\n",
    "for i in range(len(distribution)):\n",
    "    # cumulative mean\n",
    "    plt.plot(np.arange(n), np.cumsum(path == i) / np.arange(1, n + 1), label=f'proportion_{i}')\n",
    "    plt.plot(np.arange(n), distribution[i] / distribution.sum() * np.ones(n), ls='--', label=f'true_proportion_{i}')\n",
    "\n",
    "plt.xlabel('days')\n",
    "plt.title('distribution ratio per island')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1.2 MCMC algorithm\n",
    "* Pick a neighbouring island. Go with probability **P**, otherwise stay for another day.\n",
    "\n",
    "$$p = min(1, \\frac{\\text{#proposal}}{\\text{#current}})$$\n",
    "\n",
    "### Exercise\n",
    "Implement the above metropolis algorithm.\n",
    "* Loop, and in each loop\n",
    "    * Make a proposal, which is left or right\n",
    "    * Accept or reject the proposal based on the rejection rule\n",
    "    * Store all states  \n",
    "* Plot the distribution over the states\n",
    "* Plot the path of the Markov Chain that you created"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 6184 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe3f40d3438>]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 22540 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "distribution = np.array([100, 200, 30, 10, 200, 50, 60, 80, 10, 1])\n",
    "\n",
    "def generate_path(distribution, start=0):\n",
    "    current_position = start\n",
    "    while True:\n",
    "        # make a proposol (left or right)\n",
    "        proposal = (current_position + np.random.choice([-1, 1])) % len(distribution)\n",
    "        \n",
    "        # accept or reject it based on proportions\n",
    "        p = min((1, distribution[proposal] / distribution[current_position]))\n",
    "        \n",
    "        if np.random.random() < p:\n",
    "            current_position = proposal\n",
    "        yield current_position\n",
    "\n",
    "gen = generate_path(distribution, 0)\n",
    "path = np.array([next(gen) for _ in range(10000)])\n",
    "\n",
    "sns.distplot(path, kde=False, bins=len(distribution), norm_hist=True)\n",
    "plt.show()\n",
    "plt.plot(path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Things to notice\n",
    "The sample paths have become probablistic\n",
    "\n",
    "* The sample paths are different each time\n",
    "* It is independent of absolute sizes\n",
    "* We'll see much faster convergence to true distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe3f27b68d0>"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 6432 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 54524 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "distribution = np.array([50, 100, 10])\n",
    "gen = generate_path(distribution, 0)\n",
    "path = np.array([next(gen) for _ in range(1000)])\n",
    "burn_in = 200\n",
    "\n",
    "(pd.Series(path[burn_in:]).value_counts(normalize=1)\n",
    " .reindex(np.arange(len(distribution)))\n",
    " .rename('approx').to_frame()\n",
    " .assign(true=distribution / distribution.sum())\n",
    ").plot.bar()\n",
    "\n",
    "plt.figure(figsize=(16, 6))\n",
    "\n",
    "for i in range(len(distribution)):\n",
    "    # cumulative mean\n",
    "    plt.plot(np.arange(n)[burn_in:], np.cumsum(path[burn_in:] == i) / np.arange(1, n - burn_in + 1), label=f'proportion_{i}')\n",
    "    plt.plot(np.arange(n), distribution[i] / distribution.sum() * np.ones(n), ls='--', label=f'true_proportion_{i}')\n",
    "\n",
    "plt.xlabel('days')\n",
    "plt.title('distribution ratio per island')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MCMC in discrete case\n",
    "\n",
    "* Multiply the likelihood with the prior to get an unnormalized posterior\n",
    "* Sample a path with MCMC;\n",
    "    * Proposals are next or previous\n",
    "    * Accept according to relative density\n",
    "* Count visits to each state to get proportional samples of the posterior\n",
    "\n",
    "## MCMC in continuous case\n",
    "\n",
    "* Sample a chain with MCMC\n",
    "    * Get proposals via a random walk in parameter space (take steps with multivariate normal distribution)\n",
    "    * Accept or reject according to relative density of likelihood times posterior\n",
    "* These samples are proportional to the posterior distribution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Let's take a walk](http://elevanth.org/blog/2017/11/28/build-a-better-markov-chain/)\n",
    "\n",
    "*watch different mcmc implementations at the excelent blog post Markov Chains: Why Walk When You Can Flow?*"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
