{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Metropolis Hastings Algorithm\n",
    "\n",
    "*Reference -- [Youtube Video](https://www.youtube.com/watch?v=h1NOS_wxgGg)\n",
    "\n",
    "The algorithm is described below --\n",
    "\n",
    "1. Intialize $x^{(0)}$.\n",
    "2. For $i = 0$ to $N - 1$\n",
    "    1. Sample $u \\sim \\mathcal{U}(0, 1)$\n",
    "    2. Sample $x^* \\sim q(x^*|x^{(i)})$\n",
    "    3. If $u \\leq \\mathcal{A(x^{(i)}, x^*)} = \\min \\frac{p(x^*)q(x^{(i)}| x^*)}{p(x^{(i)})q(x^*|x^{(i)})}$\n",
    "       $ x^{(i+1)} = x^*$\n",
    "       else\n",
    "       $x^{(i+1)} = x^{(i)}$\n",
    "       \n",
    "# Code for Metropolis Hastings Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Import the required modules\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm, truncnorm\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEKCAYAAAD5MJl4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGKhJREFUeJzt3X+wX3V95/HnS7KBVX4UUJI2/FQIRrpdwG26LrPrZdny\no+0I21kxui2oaccRXVnZzkC0NcG1P2CVYVsLM1WEQLUpsmMBpRgoXl0dkFRF0GQh7ppIArkKEfw5\nNoH3/nHOjV/uuZf768v9kTwfM9/JuZ/v+Xy+n/Pl3vP6fj7nc76kqpAkqdeLZrsDkqS5x3CQJHUY\nDpKkDsNBktRhOEiSOgwHSVKH4aAXRJJrk7y3T20dleQHSdL+/Lkkb+1H2217dyT53X61N4nX/UCS\n7yV5bKZfe6KSPJvk5bPdD828BbPdAc0/SbYARwC7gGeAjcBNwF9Ve+NMVb19gm19G1hZVfeMtU9V\nPQocPM1uD7/eauAVVXVBT/u/0Y+2J9mPo4BLgKOq6smZfv1J8EaofZQjB01FAb9ZVYcAxwB/BlwK\nXNfvF0qyX7/bnCOOAZ6Y48EAkNnugGaH4aCpCkBV/bCqPg28AbgwyasAklyf5P3t9uFJbk/y/SRP\nJvl8W34jcDRweztt9AdJjmmnMt6aZCvwDz1lvb+vxyf5cpKnk3wqyS+0bb42yaPP6Wjy7ST/PslZ\nwHuANyT5YZKvtc/vmaZK4w+TbEmyI8kNSQ5unxvuxwVJtib5bpL3jPkGJQcnubHd79vD02xJzgDW\nA7/UHvfHRqk76nvWPndpkm+1db+R5Lye5y5M8sUkV7V1v5XkNW35d9pjuqBn/+vbKcD1bXufS3L0\nGMezMMkH22N/PMk1SfYfr7+anwwH9UVVbQC2Af92lKf/G/AocDjNdNR72joXAN8BfquqDq6qD/bU\n+XfAK4Gzhl9iRJu/C7wZWEwztfUXvd0Zo4+fBf4E+NuqOqiqThllt7cAFwCvBV4OHAR8eMQ+pwEn\nAP8BeF+SE0d7vbbeQcCxwABwQZK3VNU/AOcAj7XHPdr1k1Hfs9a3gNOq6mDgcuCvkyzqeX458ABw\nGPA3wDrgXwGvoHnfPpzkxT37v6lt53Dg68DHxzieK4DjgV9p/10CvG8C/dU8ZDionx6jOSGNtAv4\nReC4qnqmqr404vmRUxcFrK6qn1bVz8Z4rZuqalNV/RT4I+D1wxesp+lNwFVVtbWqfgKsAlb0jFoK\nWFNV/1RVD9KcTP/lyEba/d8AXFZVP6mqrcCHaE7OEzHme1ZV/6uqhtrtTwKbaQJh2Ler6sb2+s/f\nAkcCl1fVrqq6C/gnmpP7sM9U1ZeqahfwXuA1SZaM0qffB95dVU9X1Y9pphPfOF5/NT8ZDuqnJcDO\nUcr/B/B/gfXtNMelE2hr2zjP904dbQX+GfDSCfXy+f1S215v2wuA3k/mQz3bPwEOHKWdl7b1vjOi\nrdFOuqO5kjHes3Za62vtFM73gZN47rH39u+nAFX1xIiy3j7veS/bk/5OmvdhjyQvA14MfCXJziQ7\ngb+nGSnA1P4baw4zHNQXSX6V5oTyv0c+V1U/qqo/qKpXAK8DLkly+vDTYzQ53iqZo3q2j6H55PoE\n8GOak9hwv/YDXjaJdh9r2xvZ9tDou4/pibbeyLa2T6RyVf14tPesvR7wV8BFVXVoVR0KfJPpXTje\n814mOZBm9Deyn0/QBOFJVXVY+/iFdlHCeP+NNQ8ZDpqWJAcl+S2aue2bqmrjKPv8ZpJXtD/+ENhN\nc50AmpPuyHX0o53oRpb9TpJXtnPnlwOfbKdRHgEOSHJOkgXAHwILe+oNAcc+zxTU3wDvTnJse6L8\nY2BdVT37PH3raPe/GfjjJAcmOQZ4N82S33GN8Z49C7yk/feJJC9K8hbgl8drbpznfyPJv0myEPjv\nwL1V9Zx7L9r39iPA1e0ogiRLkpw5Tn81TxkOmqrbkzxNM22yCvggMNaNaScAdyf5IfAl4C+r6gvt\nc38K/FE7VXFJWzbap/sasX0TsJbmk/5C4GKAqvoBcBHNstptNCeq3imqT9KcLJ9M8o+jtP2xtu0v\n0EyT/AR41xj9GKuvw97V1v9/bXt/XVXXP8/+vUZ7zz5fVZtorl3cB+ygmVL64jhtjdfnTwBrgCeB\nU4DfGWPfS2kuht+X5CmaFVdLn6+/4x2k5q6M9z/7aZeqfYHmD3ABcEtVXZ7kUJqLXccAW4Dzq+rp\nts4qmhPFbuDiqlrflp8K3AAcANxRVf+1LV8I3Ai8mmb4+oaq6p2rlfQCSHI98GhVvW/cnbVPGXfk\n0K4WOb1d9ncycE6S5cBlwN1VdSJwD82nR9Kscz8fWEazXO+aniH8tTR3wy4FlqZZdw6wEthZVScA\nV9NcjJMkzZIJTSu1S/oA9qcZPRRwLs2wnvbf4RtxXkczR7u7qrbQLrNLshg4qF0PD81IYbhOb1u3\nAGdM6WgkTZZfj6FRTei7ldo121+huYnmL6tqQ5JFPWutdyQ5ot19CXBvT/Xtbdlunjv3u42fL+tb\nQrucrqqeSfJUksOqarRlkZL6ZIwb8KQJjxyebaeVjqQZBZzE5C7MTZbf5yJJs2hS38paVT9IMgic\nDQwNjx7aKaPvtrtt57lr0I9sy8Yq763zWLsu/eDRRg1JHAJL0hRU1aQ+dI87ckjy0iSHtNv/HPh1\nYBNwG8132wBcCNzabt9G83UDC5McR3Ob/v1VtQN4Osny9gL1BSPqXNhuv57mAveoqspHnx6rV6+e\n9T7sLQ/fS9/PufyYiomMHH4RWNted3gRzZeW3ZHkPuDmNN9muZVmhRJVtTHJzTTf8b+L5k7O4d69\ng+cuZb2zLb8OuCnJZpq11iumdDSSpL4YNxyq6iHg1FHKd9J8K+Vodf6U5uamkeVfAf7FKOU/ow0X\nSdLs8w7pfdjAwMBsd2GvsTe/l4sXH0uSGXssXnzsXv1+zhfj3iE9lySp+dRfaW/QXCKcyb+7THme\nXKNLQvX7grQkad9jOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lS\nh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUY\nDpKkDsNBktRhOEiSOsYNhyRHJrknyTeTPJTkv7Tlq5NsS/LV9nF2T51VSTYn2ZTkzJ7yU5M8mOSR\nJFf3lC9Msq6tc2+So/t9oJKkiZvIyGE3cElVnQS8Bnhnkle2z11VVae2jzsBkiwDzgeWAecA1yRJ\nu/+1wMqqWgosTXJWW74S2FlVJwBXA1f24+AkSVMzbjhU1Y6qeqDd/hGwCVjSPp1RqpwLrKuq3VW1\nBdgMLE+yGDioqja0+90InNdTZ227fQtwxhSORZLUJ5O65pDkWOBk4Mtt0TuTPJDko0kOacuWAI/2\nVNveli0BtvWUb+PnIbOnTlU9AzyV5LDJ9E2S1D8TDockB9J8qr+4HUFcA7y8qk4GdgAf6mO/RhuR\nSJJmyIKJ7JRkAU0w3FRVtwJU1fd6dvkIcHu7vR04que5I9uyscp76zyWZD/g4KraOVpf1qxZs2d7\nYGCAgYGBiRyCJO0zBgcHGRwcnFYbqarxd0puBJ6oqkt6yhZX1Y52+93Ar1bVm5K8Cvg48Gs000V3\nASdUVSW5D3gXsAH4DPDnVXVnkouAX66qi5KsAM6rqhWj9KMm0l9J/dOsJ5nJv7vg33l/JaGqJjUj\nM+7IIclpwH8GHkryNZrfkvcAb0pyMvAssAV4G0BVbUxyM7AR2AVc1HNGfwdwA3AAcMfwCifgOuCm\nJJuBJ4FOMEiSZs6ERg5zhSMHaeY5cpj/pjJy8A5pSVKH4SBJ6jAcJEkdhoMkqcNwkCR1GA6SpA7D\nQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwk\nSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1jBsOSY5Mck+S\nbyZ5KMm72vJDk6xP8nCSzyY5pKfOqiSbk2xKcmZP+alJHkzySJKre8oXJlnX1rk3ydH9PlBJ0sRN\nZOSwG7ikqk4CXgO8I8krgcuAu6vqROAeYBVAklcB5wPLgHOAa5KkbetaYGVVLQWWJjmrLV8J7Kyq\nE4CrgSv7cnSSpCkZNxyqakdVPdBu/wjYBBwJnAusbXdbC5zXbr8OWFdVu6tqC7AZWJ5kMXBQVW1o\n97uxp05vW7cAZ0znoCRJ0zOpaw5JjgVOBu4DFlXVEDQBAhzR7rYEeLSn2va2bAmwrad8W1v2nDpV\n9QzwVJLDJtM3SVL/LJjojkkOpPlUf3FV/ShJjdhl5M/TkbGeWLNmzZ7tgYEBBgYG+viykjT/DQ4O\nMjg4OK02UjX+OT3JAuDTwN9X1f9syzYBA1U11E4Zfa6qliW5DKiquqLd705gNbB1eJ+2fAXw2qp6\n+/A+VfXlJPsBj1fVEaP0oybSX0n901wynMm/u+DfeX8loarG/NA9molOK30M2DgcDK3bgDe32xcC\nt/aUr2hXIB0HHA/c3049PZ1keXuB+oIRdS5st19Pc4FbkjRLxh05JDkN+ALwEM3HhwLeA9wP3Awc\nRTMqOL+qnmrrrKJZgbSLZhpqfVv+auAG4ADgjqq6uC3fH7gJOAV4EljRXswe2RdHDtIMc+Qw/01l\n5DChaaW5wnCQZp7hMP+9kNNKkqR9iOEgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJ\nUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1\nGA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDJKlj3HBIcl2SoSQP9pStTrItyVfbx9k9z61K\nsjnJpiRn9pSfmuTBJI8kubqnfGGSdW2de5Mc3c8DlCRN3kRGDtcDZ41SflVVndo+7gRIsgw4H1gG\nnANckyTt/tcCK6tqKbA0yXCbK4GdVXUCcDVw5dQPR5LUD+OGQ1V9Efj+KE9llLJzgXVVtbuqtgCb\ngeVJFgMHVdWGdr8bgfN66qxtt28Bzph49yVJL4TpXHN4Z5IHknw0ySFt2RLg0Z59trdlS4BtPeXb\n2rLn1KmqZ4Cnkhw2jX5JkqZpwRTrXQO8v6oqyQeADwG/16c+jTYi2WPNmjV7tgcGBhgYGOjTy0rS\n3mFwcJDBwcFptZGqGn+n5Bjg9qr6led7LsllQFXVFe1zdwKrga3A56pqWVu+AnhtVb19eJ+q+nKS\n/YDHq+qIMfpRE+mvpP5pLhvO5N9d8O+8v5JQVc/7wXukiU4rhZ5P9O01hGG/DXyj3b4NWNGuQDoO\nOB64v6p2AE8nWd5eoL4AuLWnzoXt9uuBeyZzAJKk/ht3WinJJ4AB4PAk36EZCZye5GTgWWAL8DaA\nqtqY5GZgI7ALuKjno/47gBuAA4A7hlc4AdcBNyXZDDwJrOjLkUmSpmxC00pzhdNK0sxzWmn+eyGn\nlSRJ+xDDQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofh\nIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1LJjt\nDkianMWLj2VoaOtsd0N7uVTVbPdhwpLUfOqv9EJIAszk38HMv55/5/2VhKrKZOo4rSRJ6jAcJEkd\n44ZDkuuSDCV5sKfs0CTrkzyc5LNJDul5blWSzUk2JTmzp/zUJA8meSTJ1T3lC5Osa+vcm+Tofh6g\nJGnyJjJyuB44a0TZZcDdVXUicA+wCiDJq4DzgWXAOcA1aSZIAa4FVlbVUmBpkuE2VwI7q+oE4Grg\nymkcjySpD8YNh6r6IvD9EcXnAmvb7bXAee3264B1VbW7qrYAm4HlSRYDB1XVhna/G3vq9LZ1C3DG\nFI5DktRHU73mcERVDQFU1Q7giLZ8CfBoz37b27IlwLae8m1t2XPqVNUzwFNJDptivyRJfdCv+xz6\nue7seZdbrVmzZs/2wMAAAwMDfXxpSZr/BgcHGRwcnFYbUw2HoSSLqmqonTL6blu+HTiqZ78j27Kx\nynvrPJZkP+Dgqto51gv3hoMkqWvkB+fLL7980m1MdFopPPcT/W3Am9vtC4Fbe8pXtCuQjgOOB+5v\np56eTrK8vUB9wYg6F7bbr6e5wC1JmkXj3iGd5BPAAHA4MASsBv4O+CTNJ/6twPlV9VS7/yqaFUi7\ngIuran1b/mrgBuAA4I6qurgt3x+4CTgFeBJY0V7MHq0v3iGtfZ53SGuypnKHtF+fIc0zhoMmy6/P\nkCT1heEgSeowHCRJHYaDJKnD/9mPpDlmf37+lWwvvEWLjmHHji0z9nrzhauVpHlmX1it5Oqo/nK1\nkiSpLwwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaD\nJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHVMKxySbEny9SRfS3J/\nW3ZokvVJHk7y2SSH9Oy/KsnmJJuSnNlTfmqSB5M8kuTq6fRJkjR90x05PAsMVNUpVbW8LbsMuLuq\nTgTuAVYBJHkVcD6wDDgHuCZJ2jrXAiuraimwNMlZ0+yXJGkaphsOGaWNc4G17fZa4Lx2+3XAuqra\nXVVbgM3A8iSLgYOqakO73409dSRJs2C64VDAXUk2JPm9tmxRVQ0BVNUO4Ii2fAnwaE/d7W3ZEmBb\nT/m2tkySNEsWTLP+aVX1eJKXAeuTPEwTGL1G/ixJmuOmFQ5V9Xj77/eS/B2wHBhKsqiqhtopo++2\nu28HjuqpfmRbNlb5qNasWbNne2BggIGBgekcgiTtdQYHBxkcHJxWG6ma2gf7JC8GXlRVP0ryEmA9\ncDlwBrCzqq5IcilwaFVd1l6Q/jjwazTTRncBJ1RVJbkPeBewAfgM8OdVdecor1lT7a+0t2jWcczk\n38He/3p7+3klCVWV8ff8uemMHBYBn0pSbTsfr6r1Sf4RuDnJW4GtNCuUqKqNSW4GNgK7gIt6zvTv\nAG4ADgDuGC0YJEkzZ8ojh9ngyEFy5PBCvN7efl6ZysjBO6QlSR2GgySpw3CQJHUYDpKkDsNBktRh\nOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaD\nJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkjjkTDknOTvJ/\nkjyS5NLZ7o8k7cvmRDgkeRHwYeAs4CTgjUleObu92vsNDg7Odhf2Gr6X/TY42x3Y582JcACWA5ur\namtV7QLWAefOcp/2ep7Q+mPx4mM5/fTTSTIjj33D4Gx3YJ83V8JhCfBoz8/b2jJpzhsa2gqsBmqG\nHuqv/Wcs2BcvPna2D3bC5ko4zCl33HHHjP2yDD+eeOKJ2T5saR/1M2Yq2JsPEvNDqmb/k0iSfw2s\nqaqz258vA6qqrhix3+x3VpLmoaqa1JzkXAmH/YCHgTOAx4H7gTdW1aZZ7Zgk7aMWzHYHAKrqmSTv\nBNbTTHVdZzBI0uyZEyMHSdLcMi8uSCf5T0m+keSZJKeOeG5Vks1JNiU5c7b6OF8lWZ1kW5Kvto+z\nZ7tP8403cPZXki1Jvp7ka0nun+3+zDdJrksylOTBnrJDk6xP8nCSzyY5ZLx25kU4AA8B/xH4fG9h\nkmXA+cAy4Bzgmuw7C8H76aqqOrV93DnbnZlPvIHzBfEsMFBVp1TV8tnuzDx0Pc3vY6/LgLur6kTg\nHmDVeI3Mi3CoqoerajMw8sR/LrCuqnZX1RZgM80NdZocA3XqvIGz/8I8OTfNRVX1ReD7I4rPBda2\n22uB88ZrZ77/Bxh589x2vHluKt6Z5IEkH53IcFPP4Q2c/VfAXUk2JPn92e7MXuKIqhoCqKodwBHj\nVZgTq5UAktwFLOotovkleW9V3T47vdo7PN97C1wDvL+qKskHgKuAlTPfS2mP06rq8SQvowmJTe2n\nYfXPuCuR5kw4VNWvT6HaduConp+PbMvUYxLv7UcAg3hytgNH9/zs7+A0VdXj7b/fS/Ipmqk7w2F6\nhpIsqqqhJIuB745XYT5OK/XOj98GrEiyMMlxwPE0N9BpgtpflGG/DXxjtvoyT20Ajk9yTJKFwAqa\n30tNQZIXJzmw3X4JcCb+Tk5F6J4r39xuXwjcOl4Dc2bk8HySnAf8BfBS4NNJHqiqc6pqY5KbgY3A\nLuCi8saNyboyyck0K0S2AG+b3e7ML97A2XeLgE+1X5WzAPh4Va2f5T7NK0k+AQwAhyf5Ds23Qv4Z\n8MkkbwW20qzyfP52PJdKkkaaj9NKkqQXmOEgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6\n/j+pdSmcSSATCAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x9f61da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example 1 -- Normal Distribution\n",
    "\n",
    "# Initialize x_0\n",
    "x_0 = -10\n",
    "x_i = x_0\n",
    "no_samples = 50000\n",
    "samples = []\n",
    "# Target distribution p is choosen to be a Gaussian of mean 4 and scale 1\n",
    "# Proposal distribution q is choosen to be a Gaussian of mean x_1 and scale 2\n",
    "for i in range(no_samples):\n",
    "    u = np.random.rand()\n",
    "    x_star = np.random.normal(loc=x_i, scale=2)\n",
    "    if u < min(1,(norm.pdf(x_star, 4, 1)/norm.pdf(x_i, 4, 1))*(norm.pdf(x_i, x_star, 2)/norm.pdf(x_star, x_i, 2))):\n",
    "        x_i = x_star\n",
    "    samples.append(x_i)\n",
    "plt.hist(samples)\n",
    "plt.title('Distribution of samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEKCAYAAADn+anLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF95JREFUeJzt3X+QpVV95/H3ByYjUQcCmhni8MNfDD80RtCMcc2uraRA\nNAWUVeJoKajs1lZwV1Y3uzAmyrCbRHHd7JRroOKGyI9oRsQyYiQwIHZcLVCSRTFCYDAwMoM0gRH8\nGTPgd/94niF32j7Td7rv9I/h/arqmqfPPefcc+7tfj73Oc/z9KSqkCRpKvvM9wAkSQuXISFJajIk\nJElNhoQkqcmQkCQ1GRKSpCZDQntEkouS/M6I+jo0yfeSpP/+C0neNoq++/6uTvLmUfW3G8/7e0n+\nMcl9c/3cw0ry0yTPnu9xaP4sme8BaPFJcg+wHNgOPAbcBlwOfKT6G2+q6reG7Otu4MyquqFVp6ru\nBfaf5bB3PN95wHOq6vSB/l89ir53cxyHAu8CDq2qh+b6+XeDN1I9wXkkoZko4DVVdQBwOPB+4Bzg\n4lE/UZJ9R93nAnE48OACDwiAzPcANL8MCc1UAKrq+1X1l8DrgTOSHAOQ5KNJ/lu//bQkn03y3SQP\nJfnrvvwy4DDgs/1y0m8nObxf4nhbks3A5wfKBn9en5vkK0keSfLpJL/Q9/nyJPfuNNDk7iSvTHIi\n8G7g9Um+n+SW/vHHl6/S+d0k9yS5P8klSfbvH9sxjtOTbE7yQJJ3N1+gZP8kl/X17t6x/JbkeGAj\n8Ix+3n86RdspX7P+sXOS3NW3/bskpw48dkaSLyX5w77tXUle2pd/u5/T6QP1P9ovDW7s+/tCksMa\n81ma5IP93L+T5MIkT5puvFrcDAmNRFXdDGwB/vUUD/9n4F7gaXTLVO/u25wOfBv4zarav6o+ONDm\n3wBHASfueIpJfb4ZeAtwMN2S1/8eHE5jjNcCfwB8oqqWVdWxU1R7K3A68HLg2cAy4MOT6rwMOAL4\nDeC9SY6c6vn6dsuAZwJjwOlJ3lpVnwdOAu7r5z3V+ZUpX7PeXcDLqmp/4Hzgz5KsGHh8NfA14CDg\nz4ENwIuB59C9bh9O8uSB+m/s+3ka8HXgY435XAA8F3hB/+9K4L1DjFeLmCGhUbqPbsc02Xbgl4Bn\nVdVjVfXlSY9PXtIo4Lyq+nFV/aTxXJdX1e1V9WPgPcDrdpzYnqU3An9YVZur6kfAWmDNwFFMAeuq\n6p+r6la6neqvTO6kr/964Nyq+lFVbQb+J91OehjN16yqPlVVE/32J4FNdMGww91VdVl/fugTwCHA\n+VW1vaquA/6Zbie/w+eq6stVtR34HeClSVZOMaZ/B7yzqh6pqh/SLTO+YbrxanEzJDRKK4FtU5T/\nD+BbwMZ++eOcIfraMs3jg0tKm4GfA54+1Ch37Rl9f4N9LwEGP6lPDGz/CHjqFP08vW/37Ul9TbXz\nncoHaLxm/XLXLf3SzneB57Hz3AfH92OAqnpwUtngmB9/Lfud/za61+FxSX4ReDLwt0m2JdkG/BXd\nkQPM7D3WImBIaCSS/CrdjuX/Tn6sqn5QVb9dVc8BTgbeleQVOx5udDndVTWHDmwfTvdJ9kHgh3Q7\nsx3j2hf4xd3o976+v8l9T0xdvenBvt3kvrYO07iqfjjVa9afL/gIcFZVHVhVBwLfZHYnmB9/LZM8\nle5ocPI4H6QLxOdV1UH91y/0Fy9M9x5rETMkNCtJliX5Tbq178ur6rYp6rwmyXP6b78PPEp3HgG6\nne/k6/Cn2uFNLntTkqP6tfXzgU/2yyt3AvslOSnJEuB3gaUD7SaAZ+5iaerPgXcmeWa/w/x9YENV\n/XQXY/sZff0rgN9P8tQkhwPvpLtUeFqN1+ynwFP6fx9Msk+StwLPn667aR5/dZJ/lWQp8N+BG6tq\np3s3+tf2/wDr+6MKkqxMcsI049UiZ0hopj6b5BG65ZS1wAeB1g1uRwDXJ/k+8GXgj6rqi/1j7wPe\n0y9hvKsvm+rTfk3avhy4lO6T/1LgbICq+h5wFt3luFvodliDS1efpNtpPpTkb6bo+0/7vr9It3zy\nI+AdjXG0xrrDO/r2/9D392dV9dFd1B801Wv211V1O925jZuA++mWmr40TV/TjfnjwDrgIeBY4E2N\nuufQnTS/KcnDdFdordrVeKebpBa+DPOfDqW7eeoRuk8G26tqdZID6U6KHQ7cA5xWVY/09dfS7TAe\nBc6uqo19+XHAJcB+wNVV9Z/68qXAZcCL6A5rX19Vg2u5kvaAJB8F7q2q905bWU9Iwx5J/BQYq6pj\nq2rHVRTnAtdX1ZHADXSfJkl3nfxpwNF0l/ldOHBofxHd3bWrgFXprlsHOBPYVlVHAOvpTtpJkubZ\nsCGRKeqeQne4T//vjht6TqZbw320qu6hvzwvycHAsv56euiOHE6doq8rgeN3ZxKSZsw/u6FdGvZv\nNxVwXZLHgD+uqj8BVgxcq31/kuV93ZXAjQNtt/Zlj7Lz2vAW/uVywJX0l+FV1WNJHk5yUFVNdTml\npBFp3MgnPW7YkHhZVX2nv6phY5I72L0TeLvLvxcjSQvAUCFRVd/p//3HJH9Bd3fnRJIVVTXRLyU9\n0Fffys7XsB/Sl7XKB9vc11/Xvv9URxFJPDSWpBmoqhl9+J72nESSJ/fXi5PkKcAJwDeAq+j+dg7A\nGcBn+u2r6P6MwdIkz6K7/f+rVXU/8EiS1f2J7NMntTmj334d3YnwKVXVXvt13nnnzfsYnN/Mvvqf\nzjn8mtvfhb35vXsizG82hjmSWAF8uv8UvwT4WFVt7K8xvyLdX8/cTHdFE1V1W5Ir6P6Pge10d4bu\nGOXb2fkS2Gv68ouBy5NsortWe82sZiVJGolpQ6Kq7gZeOEX5Nrq/gjlVm/fR3SQ1ufxvgV+eovwn\n9CEjSVo4vON6ARkbG5vvIexRe/v89mZ7+3u3t89vNoa643qhSFKLabx64uhOs83lz2ZmvdasJ44k\n1J46cS1JeuIyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQkSU2GhCSpyZCQ\nJDUZEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLUZEhIkpoMCUlS\nkyEhSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQkSU1Dh0SSfZL8vyRX\n9d8fmGRjkjuSXJvkgIG6a5NsSnJ7khMGyo9LcmuSO5OsHyhfmmRD3+bGJIeNaoKSpJnbnSOJs4Hb\nBr4/F7i+qo4EbgDWAiQ5BjgNOBo4CbgwSfo2FwFnVtUqYFWSE/vyM4FtVXUEsB74wAznI0kaoaFC\nIskhwKuBPxkoPgW4tN++FDi13z4Z2FBVj1bVPcAmYHWSg4FlVXVzX++ygTaDfV0JHL/7U5Ekjdqw\nRxL/C/gvQA2UraiqCYCquh9Y3pevBO4dqLe1L1sJbBko39KX7dSmqh4DHk5y0PDTkCTtCUumq5Dk\nNcBEVX0tydguqtYuHttdaT2wbt26x7fHxsYYGxsb4dNK0uI3Pj7O+Pj4SPpK1a737Un+AHgT8Cjw\n88Ay4NPAi4Gxqprol5K+UFVHJzkXqKq6oG9/DXAesHlHnb58DfDyqvqtHXWq6itJ9gW+U1XLJw2F\nJDXdeKX50J12m8ufzeDvgoaVhKpqfvjelWmXm6rq3VV1WFU9G1gD3FBVbwY+C7ylr3YG8Jl++ypg\nTX/F0rOA5wJf7ZekHkmyuj+RffqkNmf026+jOxEuSZpn0y437cL7gSuSvI3uKOE0gKq6LckVdFdC\nbQfOGvj4/3bgEmA/4OqquqYvvxi4PMkm4CG6MJIkzbNpl5sWEpebtFC53KSFbI8uN0mSnrgMCUlS\nkyEhSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQkSU2GhCSpyZCQJDUZ\nEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLUZEhIkpoMCUlSkyEh\nSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqmjYkkjwpyVeS3JLkG0nO68sP\nTLIxyR1Jrk1ywECbtUk2Jbk9yQkD5ccluTXJnUnWD5QvTbKhb3NjksNGPVFJ0u6bNiSq6ifAK6rq\nWOCFwElJVgPnAtdX1ZHADcBagCTHAKcBRwMnARcmSd/dRcCZVbUKWJXkxL78TGBbVR0BrAc+MKoJ\nSpJmbqjlpqr6Ub/5JGAJUMApwKV9+aXAqf32ycCGqnq0qu4BNgGrkxwMLKuqm/t6lw20GezrSuD4\nGc1GkjRSQ4VEkn2S3ALcD1zX7+hXVNUEQFXdDyzvq68E7h1ovrUvWwlsGSjf0pft1KaqHgMeTnLQ\njGYkSRqZYY8kftovNx1Cd1TwPLqjiZ2qjXBcmb6KJGlPW7I7lavqe0nGgVcBE0lWVNVEv5T0QF9t\nK3DoQLND+rJW+WCb+5LsC+xfVdumGsO6dese3x4bG2NsbGx3piBJe73x8XHGx8dH0leqdn0AkOTp\nwPaqeiTJzwPXAu8HXk53svmCJOcAB1bVuf2J648BL6FbRroOOKKqKslNwDuAm4HPAR+qqmuSnAU8\nv6rOSrIGOLWq1kwxlppuvNJ86K7NmMufzeDvgoaVhKqa0QrNMEcSvwRcmmQfuuWpT1TV1f0O/4ok\nbwM2013RRFXdluQK4DZgO3DWwJ797cAlwH7A1VV1TV9+MXB5kk3AQ8DPBIQkae5NeySxkHgkoYXK\nIwktZLM5kvCOa0lSkyEhSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQk\nSU2GhCSpyZCQJDUZEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLU\nZEhIkpoMCUlSkyEhSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQkSU2G\nhCSpadqQSHJIkhuSfDPJN5K8oy8/MMnGJHckuTbJAQNt1ibZlOT2JCcMlB+X5NYkdyZZP1C+NMmG\nvs2NSQ4b9UQlSbtvmCOJR4F3VdXzgJcCb09yFHAucH1VHQncAKwFSHIMcBpwNHAScGGS9H1dBJxZ\nVauAVUlO7MvPBLZV1RHAeuADI5mdJGlWpg2Jqrq/qr7Wb/8AuB04BDgFuLSvdilwar99MrChqh6t\nqnuATcDqJAcDy6rq5r7eZQNtBvu6Ejh+NpOSJI3Gbp2TSPJM4IXATcCKqpqALkiA5X21lcC9A822\n9mUrgS0D5Vv6sp3aVNVjwMNJDtqdsUmSRm/okEjyVLpP+Wf3RxQ1qcrk72cj01eRJO1pS4aplGQJ\nXUBcXlWf6Ysnkqyoqol+KemBvnwrcOhA80P6slb5YJv7kuwL7F9V26Yay7p16x7fHhsbY2xsbJgp\nSNITxvj4OOPj4yPpK1XTHwAkuQx4sKreNVB2Ad3J5guSnAMcWFXn9ieuPwa8hG4Z6TrgiKqqJDcB\n7wBuBj4HfKiqrklyFvD8qjoryRrg1KpaM8U4apjxSnOtuzZjLn82g78LGlYSqmpGKzTThkSSlwFf\nBL5B91tQwLuBrwJX0B0BbAZOq6qH+zZr6a5Y2k63PLWxL38RcAmwH3B1VZ3dlz8JuBw4FngIWNOf\n9J48FkNCC5IhoYVsj4bEQmJIaKEyJLSQzSYkvONaktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQ\nkCQ1GRKSpCZDQpLUZEhIkpoMCUlSkyEhSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJ\nUpMhIUlqMiQkSU2GhCSpyZCQJDUZEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1\nGRKSpCZDQpLUZEhIkpoMCUlSkyEhSWqaNiSSXJxkIsmtA2UHJtmY5I4k1yY5YOCxtUk2Jbk9yQkD\n5ccluTXJnUnWD5QvTbKhb3NjksNGOUFJ0swNcyTxUeDESWXnAtdX1ZHADcBagCTHAKcBRwMnARcm\nSd/mIuDMqloFrEqyo88zgW1VdQSwHvjALOYjSRqhaUOiqr4EfHdS8SnApf32pcCp/fbJwIaqerSq\n7gE2AauTHAwsq6qb+3qXDbQZ7OtK4PgZzEOStAfM9JzE8qqaAKiq+4HlfflK4N6Belv7spXAloHy\nLX3ZTm2q6jHg4SQHzXBckqQRWjKifmpE/QBkVw+uW7fu8e2xsTHGxsZG+NSStPiNj48zPj4+kr5m\nGhITSVZU1US/lPRAX74VOHSg3iF9Wat8sM19SfYF9q+qba0nHgwJSdLPmvwB+vzzz59xX8MuN4Wd\nP+FfBbyl3z4D+MxA+Zr+iqVnAc8FvtovST2SZHV/Ivv0SW3O6LdfR3ciXJK0AKRq1ytFST4OjAFP\nAyaA84C/AD5JdwSwGTitqh7u66+lu2JpO3B2VW3sy18EXALsB1xdVWf35U8CLgeOBR4C1vQnvaca\nS003Xmk+dJ995vJnM/i7oGEloap2uZTfbLuYftAMCS1UhoQWstmEhHdcS5KaDAlJUpMhIUlqMiQk\nSU2GhCSpyZCQJDUZEpKkJkNCktRkSEiSmgwJSVKTISFJajIkJElNhoQkqcmQkCQ1GRKSpCZDQpLU\nZEhIkpoMCUlSkyEhSWoyJCRJTYaEJKnJkJAkNRkSkqQmQ0KS1GRISJKaDAlJUpMhIUlqMiQkSU2G\nhCSpyZCQJDUZEpKkJkNCktRkSEiSmgwJSVKTISFJalowIZHkVUn+PsmdSc6Z7/FIkhZISCTZB/gw\ncCLwPOANSY6a31HNvfHx8fkewh61t89vb7a3v3d7+/xmY0GEBLAa2FRVm6tqO7ABOGWexzTn9vYf\n1L19fnuzvf2929vnNxsLJSRWAvcOfL+lL5MkzaMl8z2AhWr79u289rWvnbPnW758OYceeuicPZ8k\nDSNVNd9jIMmvAeuq6lX99+cCVVUXTKo3/4OVpEWoqjKTdgslJPYF7gCOB74DfBV4Q1XdPq8Dk6Qn\nuAWx3FRVjyX5D8BGuvMkFxsQkjT/FsSRhCRpYVooVzftZNgb65L8apLtSebuDPMIDDO/JGNJbkny\nd0m+MNdjnKnp5pZk/yRXJflakm8kecs8DHPGklycZCLJrbuo86Ekm/o5vnAuxzcb080tyRuTfL3/\n+lKSX57rMc7GMO9dX2+x7leG+dnc/f1KVS2oL7rgugs4HPg54GvAUY16nwf+EnjtfI97lPMDDgC+\nCazsv3/6fI97hHNbC7xvx7yAh4Al8z323ZjjrwMvBG5tPH4S8Ll++yXATfM95hHO7deAA/rtVy2m\nuQ0zv77OotyvDPn+zWi/shCPJIa9se4/AlcCD8zl4EZgmPm9EfhUVW0FqKoH53iMMzXM3ApY1m8v\nAx6qqkfncIyzUlVfAr67iyqnAJf1db8CHJBkxVyMbbamm1tV3VRVj/Tf3sQiu5dpiPcOFu9+ZZj5\nzWi/shBDYtob65I8Azi1qi4CZnRZ1zwa5sbBVcBBSb6Q5OYkb56z0c3OMHP7MHBMkvuArwNnz9HY\n5srk12Ari2xnOqR/C/zVfA9ilBb5fmUYM9qvLIirm2ZgPTC43r23vaFLgOOAVwJPAW5McmNV3TW/\nwxqJE4FbquqVSZ4DXJfkBVX1g/kemIaT5BXAW+mWN/Ym7lcajRaarcBhA98f0pcNejGwIUno1rVP\nSrK9qq6aozHOxjDz2wI8WFX/BPxTki8Cv0K33r+QDTO3twLvA6iqbyW5GzgK+Js5GeGetxUYvHV+\nqtdg0UryAuAjwKuqarqlm8VmMe9XhjGj/cpCXG66GXhuksOTLAXWADu9SVX17P7rWXTrh2ctojdy\n2vkBnwF+Pcm+SZ5MdwJ0Mdw3MszcNgO/AdCv1a8C/mFORzl7of0p8yrgdHj8Lwk8XFUTczWwEWjO\nLclhwKeAN1fVt+Z0VKPTnN8i36/ssKufzRntVxbckUQ1bqxL8u+7h+sjk5vM+SBnYZj5VdXfJ7kW\nuBV4DPhIVd02j8MeypDv3e8Blwxcpvdfq2rbPA15tyX5ODAGPC3Jt4HzgKX8y3t3dZJXJ7kL+CHd\nkdOiMN3cgPcABwEX9p+2t1fV6vka7+4aYn6DFtV+BYb62ZzRfsWb6SRJTQtxuUmStEAYEpKkJkNC\nktRkSEiSmgwJSVKTISFJajIkJElNhoQkqen/A0iXqD162ev9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa24ab38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example 2 -- Truncated Normal Distribution\n",
    "\n",
    "# Initialize x_0\n",
    "x_0 = 1\n",
    "x_i = x_0\n",
    "no_samples = 50000\n",
    "samples = []\n",
    "# Target distribution p is choosen to be a Truncated Gaussian of a = 4 b = 6 mean 4 and scale 5\n",
    "# Proposal distribution q is choosen to be a Gaussian of mean x_1 and scale 2\n",
    "for i in range(no_samples):\n",
    "    u = np.random.rand()\n",
    "    x_star = np.random.normal(loc=x_i, scale=3)\n",
    "    if u < min(1,(truncnorm.pdf(x_star, 4, 6, 5, 5)/(truncnorm.pdf(x_i, 4, 6, 5, 5)+.001))*(norm.pdf(x_i, x_star, 3)/norm.pdf(x_star, x_i, 3))):\n",
    "        x_i = x_star\n",
    "    samples.append(x_i)\n",
    "plt.hist(samples)\n",
    "plt.title('Distribution of samples')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}