{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PDF Sampling: MCMC Metropolis-Hastings\n",
    "\n",
    "We use a proposal distribution $Q(x' | x^{(t)})$ to generate a proposed next state $x'$ based on current state $x^{(t)}$.\n",
    "\n",
    "<strong>Metropolis-Hastings Algorithm:</strong>\n",
    "0. Initiate state $x^{(t=0)}$\n",
    "1. Generate a proposed state from $Q(x' | x^{(t)})$\n",
    "2. Compute $a = \\frac{f^{*}(x')}{f^{*}(x^{(t)})} \\frac{Q(x^{(t)} | x')}{Q(x' | x^{(t)})}$\n",
    "3. If $a \\ge 1$ accept proposed state $x'$\n",
    "4. Otherwise, generate $u$ from a uniform distribution in interval $[0, 1]$ and accept the proposed state if $a \\ge u$.\n",
    "5. If the state is accepted, we set $x^{(t+1)} = x'$\n",
    "6. If the state is rejected, stick to the previous state and set $x^{(t+1)} = x^{(t)}$\n",
    "\n",
    "Note that if $Q(x' | x^{(t)})$ is symmetric, the ratio of $\\frac{Q(x^{(t)} | x')}{Q(x' | x^{(t)})}$ is always 1.\n",
    "\n",
    "> * SEAYAC Workshop, *MC methods in astronomy*, Tri L. Astraatmadja (MPIA Heidelberg)\n",
    "  Krabi, 4 December 2015"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "from scipy.integrate import quad\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pdftarget(x, norm=1):\n",
    "    return np.exp(0.4*(x-0.4)*(x-0.4) - 0.08*x*x*x*x)/norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sampling_mcmc_mh(xt, stepsize, nsamp):\n",
    "    samples = np.empty(nsamp)\n",
    "    accept = np.empty(nsamp)\n",
    "    for i in range(nsamp):\n",
    "        xprime = xt + stepsize*np.random.normal() # gaussian proposal distribution\n",
    "        a = pdftarget(xprime)/pdftarget(xt) # symmetric -> gaussian\n",
    "        if a >= 1.0:\n",
    "            xt = xprime\n",
    "            accept[i] = 1\n",
    "        else:\n",
    "            u = np.random.random()\n",
    "            if a >= u:\n",
    "                xt = xprime\n",
    "                accept[i] = 1\n",
    "            else:\n",
    "                accept[i] = 0 # reject xprime, xt = xt\n",
    "        \n",
    "        samples[i] = xt\n",
    "        \n",
    "    return samples, accept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's try it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 0\n",
    "sigma = 2\n",
    "nsamp = 10000\n",
    "\n",
    "samples, accept = sampling_mcmc_mh(mu, sigma, nsamp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9xvHPNwlb2TchsoU1EFD2TdSkLCJghVYthKKI\nFexCxXpr7dXqRXrr0sV6rQsu1F0odcNakEVIgSoQdhPQsMi+SkA2DYH87h9JSIgJGcIkZ87M8369\nzss5MydnHgZ4/PGbs5hzDhER8acorwOIiEjZqcRFRHxMJS4i4mMqcRERH1OJi4j4mEpcRMTHAipx\nM7vWzD4zswwzu7eY1xPN7IiZrc5bfhv8qCIiUlRMaRuYWRTwFDAA2AOkmtks59xnRTZd7Jy7vhwy\niohICQIZifcCNjnntjvnsoEZwPBitrOgJhMRkVIFUuJNgJ2F1nflPVdUXzNba2b/MrOEoKQTEZHz\nKnU6JUCrgObOuZNmNgR4D2gXpH2LiEgJAinx3UDzQutN8547yzl3vNDjOWb2jJnVc85lFt7OzHSh\nFhGRMnDOFTtlHch0SirQxsxamFllYBTwfuENzKxRoce9ACta4IWCeL78z//8j+cZQmXRZ6HPQp9F\n6H8W51PqSNw5d8bMJgLzyC39ac65jWZ2R+7L7nngRjP7KZANfA2MLG2/IiJy8QKaE3fOfQjEF3nu\nuUKPnwaeDm40EREpTUSesZmUlOR1hJChz6KAPosC+iwKhPpnYaXNtwT1zcxcRb6fiEg4MDPcRXyx\nKSIiIUolLiLiYypxEREfU4mLiPiYSlxExMdU4iIiPqYSFxHxMZW4iIiPqcRFRHxMJS4i4mMqcRER\nH1OJi4j4mEpcRMTHgnWPTQljhw4VPDaDevW8yyIi51KJS6mmToXs7NzHVarAf/+3t3lEpICmU0RE\nfEwlLiLiYypxEREfU4mLiPiYSlxExMdU4iIiPqZDDOWinDwJqakF67VrQ5cu3uURiTQqcbkoJ0/C\nokUF63FxKnGRiqTpFBERH1OJi4j4mEpcRMTHVOIiIj6mEhcR8TGVuIiIj6nERUR8TCUuIuJjKnER\nER9TiYuI+JhKXETExwIqcTO71sw+M7MMM7v3PNv1NLNsM/tB8CKKiEhJSi1xM4sCngIGAx2BZDNr\nX8J2jwJzgx1SRESKF8hIvBewyTm33TmXDcwAhhez3S+At4ADQcwnIiLnEUiJNwF2FlrflffcWWZ2\nKTDCOfcsYMGLJyIi5xOs64k/ARSeKy+xyCdPnnz2cVJSEklJSUGKICISHlJSUkhJSQlo20BKfDfQ\nvNB607znCusBzDAzAxoAQ8ws2zn3ftGdFS5xERH5tqID3IceeqjEbQMp8VSgjZm1APYCo4Dkwhs4\n51rlPzazl4B/FlfgIiISXKWWuHPujJlNBOaRO4c+zTm30czuyH3ZPV/0R8ohp4iIFCOgOXHn3IdA\nfJHnnith29uCkEtERAKgMzZFRHxMJS4i4mPBOsRQwkxWVhaffvoplSpVIienExDtdSQRKYZKXL7l\ntdde49e//jUNGzbk1KlT7N//DYMG/ZmEhBu8jiYiRajE5RyPPvooL774Ih988AHdu3cHYPz4pbz9\n9lgOHkznmmse9DihiBSmEpezZs2axTPPPMPy5cuJjY09+3xc3JXcdtt/eOWV71K9egNmzPjZ2dey\nsrxIKiL5VOICQGZmJhMmTOC99947p8Dz1ajRmOTkfzJt2hXExvagSZNeHqQUkaJ0dIoA8OCDD3LD\nDTfQt2/fErepV68Nw4Y9w7vv3szp0xqCi4QClbiwa9cu3nzzTaZMmVLqtgkJN9KgQQf+85/HKiCZ\niJRGJS48/vjjjBs3jgYNGgS0/ZAhT7J8+f/x1Vc7S99YRMqVSjzCHT58mJdffpm777474J+pXbs5\n3bqNZ/Hi35VjMhEJhEo8wk2fPp3BgwfTpEmT0jcupF+/X7Nx4ztkZm4up2QiEgiVeIT729/+xm23\nXfg1y6pVq0fv3nfy73+XPo8uIuVHJR7B1q1bx8GDB+nfv3+Zfr537zvJyPiAo0d3BTmZiARKJR7B\n3njjDcaMGUN0dNmui1K1ah06dx7L8uVPBjmZiARKJR6hnHO888473HjjjRe1nz59JrFmzTSyso4G\nKZmIXAiVeIRav349Z86coUuXLhe1nzp14mjVahBr1rwUpGQiciFU4hHqnXfe4YYbbiD33tYXp2fP\nn7Fq1XM4pzvziVQ0lXiEmjVrFiNGjAjKvpo3vwqAHTuWBGV/IhI4lXgE2rdvH9u3b6dPnz5B2Z+Z\n0aPHT1i58tmg7E9EAqerGEagBQsW0L9/f2Jiiv/t37ABTp0qWM/JKX2fl19+M4sWPcjevQdYu/aS\ns8/XrQstWlxsYhEpiUo8As2bN49rrrmmxNfnz4fDhy9sn9Wq1aV9+xEsWvQqWVm/Ovt8p04qcZHy\npOmUCOOcK7XEy6pLl1tZt+4VfcEpUoFU4hHm008/pWbNmrRs2TLo+27R4mqyso6xb9+aoO9bRIqn\nEo8wCxYsYNCgQeWyb7MoOne+hbVrXymX/YvIt6nEI8zixYtJTEwst/137nwLaWnTOXPmVOkbi8hF\nU4lHkJycHJYuXcqVV15Zbu9Rr14b6tdvx6ZNc8rtPUSkgEo8gnz++efUqlXrgq8dfqE6dx7LunWa\nUhGpCDrEMIIsWbKkXEfh+Tp2/CHz5/+Kr7/OBOqd89rGjbB6dcF6166QkFDukUTClkbiEWTJkiVc\nddVV5f4+VavWpnXra9iw4e1vvXbkCGzaVLBc6PHoInIulXgEKe/58MI6dUomLW16hbyXSCRTiUeI\n3bt3c+zYMdq3b18h79e27VD27VtDZuaeCnk/kUilOfEIsWLFCnr37l3spWf//nfYU6hrjx27+PeL\nialKfPxwPvlkJhMm3HXxOxSRYmkkHiFWrFhBz549i33txAn46quCJZALXgWiU6dkli7VlIpIeVKJ\nR4jU1FR69epVoe/ZqtUADh7cxpYtWyr0fUUiSUAlbmbXmtlnZpZhZvcW8/r1ZrbOzNaY2Qoz6xf8\nqFJWOTk5rFy5ssSReHmJioqhT58bmT5do3GR8lJqiZtZFPAUMBjoCCSbWdFvxxY45zo757oCPwZe\nDHpSKbNNmzZRt25dGjZsWOHvfeWVo5k+fbqubChSTgIZifcCNjnntjvnsoEZwPDCGzjnThZarQEE\naVZVgiE1NbXCR+H52rXry/Hjx/n00089eX+RcBdIiTcBdhZa35X33DnMbISZbQT+CdwWnHgSDCtW\nrKjw+fB8UVFRjBo1SlMqIuUkaF9sOufec851AEYA/xus/crF83IkDpCcnMyMGTM0pSJSDgI5Tnw3\n0LzQetO854rlnFtqZq3MrJ5zLrPo65MnTz77OCkpiaSkpIDDyoU7deoU69evp3v37p5l6Ny5M1Wr\nVmXZsmVAX89yiPhFSkoKKSkpAW0bSImnAm3MrAWwFxgFJBfewMxaO+e25D3uBlQursDh3BKX8rdh\nwwbi4uKoUaOGZxnMjOTkZKZPn05yskpcpDRFB7gPPfRQiduWOp3inDsDTATmAenADOfcRjO7w8wm\n5G12g5mlmdlq4K/AD8seX4Jp7dq1dO3a1esYJCcnM3PmTE6fPu11FJGwEtBp9865D4H4Is89V+jx\nH4A/BDeaBMPatWvp0qWL1zFo27YtzZo1Y/XqRUD53B5OJBLpjM0wFyolDjB69GjmzXvT6xgiYUUl\nHsacc6xdu5bOnTt7HQWAkSNHsmTJLE6f/sbrKCJhQyUexrZv30716tU9OVOzOJdeeilt23Zh06bZ\nXkcRCRu6FG0YC4WplL17YU6heyZfdtloUlOn06HDD7wLJRJGVOJhLBRK/NCh3CVfo0Y3sGXLf5GV\ndZQqVWp5F0wkTGg6JYyFQokXVa1aXeLikti48V2vo4iEBZV4GAvFEgfdf1MkmFTiYerw4cMcOnSI\n1q1bex3lW9q1+x67di3jxIkDXkcR8T2VeJhat24dl19+OVFRofdbXLlyddq1G0Z6+j+8jiLie6H3\nN1yCIpSODy9Op06jSUvTiT8iF0slHqbS0tK4/PLLvY5RotatB/Hll5+zd+82j5OI+JtKPEylpaXR\nqVMnr2OUKDq6MgkJN7JgwQyvo1SYP/2pYHnmGa/TSLhQiYehnJwc0tPT6dixo9dRzqtTp2QWLIic\no1SOHy9YTp4sfXuRQKjEw9COHTuoXbs2devW9TrKebVocRVHjx4iPT3d6ygivqUSD0OhPpWSzyyK\ngQOTef31172OIuJbKvEw5IeplHxDhozltdde48yZM15HEfEllXgY8stIHKBVq07Exsby0UcfeR1F\nxJdU4mHITyUOcOutt/Lyyy97HUPEl3QVwzBz+vRpPv/8c6pXT2Dbttzn6tWDWiF8wcBRo0Zx//33\nc+TIEbKz63DiRMFrjRpBtWreZRMJdSrxMLNlyxbq1GnMzJnVzz43ZAj07u1hqFLUr1+fQYMGMXPm\nTBo0mMD69QWv3XILtGrlXTaRUKfplDCTnp5Oy5b+mUrJpykVkbJRiYeZtLQ0WrXyX4kPHjyYL774\ngt27P/c6ioivaDolzKSlpdG+/YjzbrN5M2RmFqwfO1bOoc5j505YsQIghgEDxjB79iv07v2wd4FE\nfEYlHmbS0tIYMuR+tm8veZvVq2HDhorLdD6ffZa7ANSuPY5ZswbSo8dDREdX8jaYiE9oOiWMZGVl\nsXXrVpo1i/c6Spk0bJhAvXqtycj4wOsoIr6hEg8jGRkZtGzZkipVqnodpcy6d/8Jq1ZN9TqGiG9o\nOiWMlHSST1oaHDxYsL53bwWGukAJCTcwd+5dZGZuoV691ixbdu7UT2Ii1KzpXb6KcuBA/ncFuZo1\ngxC+x4d4SCUeRkoq8Z07cxc/iImpSufOY1m9+gUGDnyUjIxzX+/VKzJK/KuvYOXKgvUzZ1TiUjxN\np4SRtLQ031z46ny6d5/A2rUvcfp0ltdRREKeSjyM+O2aKSWpX78dl1zSic8+e9frKCIhTyUeJk6c\nOMGePXto06aN11GConv3n7By5bNexxAJeSrxMLFx40bi4+OJiQmPrznatx/B4cNb2bt3jddRREKa\nSjxMhMt8eL7o6Er07DmRZcv+4nUUkZCmEg8T4TIfXlj37hPIyPgnx46F8DGRIh5TiYeJcCzxatXq\n0qnTaFJTn/Y6ikjIUomHCT/dV/NC9O59J6tWPU929tdeRxEJSQGVuJlda2afmVmGmd1bzOujzWxd\n3rLUzC4LflQpyZEjRzh8+DBxcXFeRwm6Bg3iadKkF+vWvep1FJGQVGqJm1kU8BQwGOgIJJtZ+yKb\nbQWuds51Bv4XeCHYQaVkGzZsoEOHDkRFhec/rPr1u5ePP/4DOTmnvY4iEnIC+VvfC9jknNvunMsG\nZgDDC2/gnFvmnPsqb3UZ0CS4MeV8wnUqJV+LFldRq1ZT0tJmeB1FJOQEUuJNgMJX3tjF+Uv6dmDO\nxYSSCxPuJQ5w1VX3s2TJw+Tk5HgdRSSkBPXf32b2XWAc8K15cyk/6enpYXdkSlGtWg2icuUazJ6t\nU/FFCgvk9L7dQPNC603znjuHmV0OPA9c65w7XNLOJk+efPZxUlISSUlJAUaVkoTbiT7FMTOuvvq3\nPPHEZMaN+wFm9q1tnIOtWwvWY2KgRYuS93n8OOzfX7BeqxY0bBjE0CJllJKSQkpKSkDbBlLiqUAb\nM2sB7AVGAcmFNzCz5sDbwM3OuS3n21nhEpeLl5mZyYkTJ2jWrJnXUcpdu3bf49NPH2TWrFmMGPHt\n+4ieOQOvvVawXqcO3HVXyfvbuRP+/veC9Z49YdiwIAYWKaOiA9yHHnqoxG1LnU5xzp0BJgLzgHRg\nhnNuo5ndYWYT8jZ7AKgHPGNma8xsRQm7kyBLT08nISGh2JFpuDEz7rvvYe6//37OnDnjdRyRkBDQ\nnLhz7kPnXLxzrq1z7tG8555zzj2f93i8c66+c66bc66rc65XeYaWApHwpWZhAwYMoX79+rz6qo4b\nFwHd2cf3IuFLzcLMjMcee4yRI0eSnJxM1ar+vZ+ol1avhqxC99zo3RvC9DSDsKffNp+LhC81i+rb\nty/dunXj6ad1TZWyWrwY5s4tWHTkpn+pxH0u0qZT8j388MM89thjfPnll15HEfGUStzHDh48yKlT\np7j00ku9jlLhEhISGD16NPfff7/XUUQ8pRL3sfz58Eg4MqU4kydP5v3332dl4dvCi0QYlbiPRepU\nSr46derw8MMPM3HiRJ2OLxFLJe5jkV7iAGPHjiUqKoqpU6d6HUXEEypxH4vEI1OKioqK4sUXX+TB\nBx9k+/btXscRqXA6TtynnHMaiedJSEjg7rvv5ic/Gc8VV8wt83cEGzbAnj0F6927Q7duQQopUk40\nEvep/fv3Y2Y0atTI6ygh4Z577uHQoUOsXv1imfdx4gTs3l2wHDsWxIAi5UQl7lP5o/BIPTKlqEqV\nKvHSS6+ycOF9HDy40es4IhVGJe5Tmkr5to4dO9K//8O89dZI3VhZIobmxH1g/354sdAsQatWuV9q\ndu7c2btQIapbt9v54ouPmDv3l4wZoyNWyuLrr+HxxwvWGzSAO+7wLo+cn0biPuAcZGcXLGfORN6F\nrwJlZlx33XNs25bCxx8/73Uc3yr85y072+s0cj4qcR9yzpGWlqYSL0HVqrVJTv4nc+Y8EPDdUUT8\nStMpPrRv33Zq1KhJdHR9jhz59utfazqY+vXbcsst0xk1ahSLFy+mXbt2XkcKqm++yV3yVasGVaoU\nrJ84ce4IukaN3NvVSfjRb6sPrV69nu9853KeeMLrJKGtbdv+PPLIIwwaNIjFixfT4nw33PSZFStg\n4cKC9aFDoVehW7F88AFsLHSQzo9/DBFwB7+IpBL3of3719Oo0eVex/CFcePGcfToUQYOHMjixYuJ\njY31OpJIUGlO3If2719Ho0Y6MiVQkyZNYty4cSQmJrJt2zav44gElUbiPrR//3oSEyd7HSPkZWdD\nRkbu4xtvvI+TJ2vQt+9VvPDCHNq163TOKfZ+d+BAwa8V4Phx77JIxVKJ+0x29km++mon9euH1xd1\n5eHECXjzzYL1mJg7ueKKBowc2Z8RI16mbduh3oULspUrcxeJPCpxnzlwII0GDdoTHV3J6yi+dNll\no6lduwVvvTWSbt3Gk5j4AGb+mFU8ffo0e/bsYf/+/aSlZbNjRw7VqtWjZs0mVK1a2+t44hGVuM/o\nS82L17x5P8aPT+Xtt0fxxRcfcf3106hfv63Xsc5x6tQptm37hN27l7NnTypTp64mM3Mnl1xyCY0b\nN+abbypz9Khx8uQhjh7dRfXqDWnatC/x8dfTrt33qFy5ute/BKkgKnGfUYkHR82asdxyy0JWrPgr\n06b15Yor7qFPn0nExFQ9u8327bB0acHPtGkDjRsXrKenw+HDBetduuQej10Wzjk2b97M3LlzmTt3\nLikpi6lVK55mza4gPn44Awb8jvHjW1G5cmUg9271+YcYOpfDoUMZ7NixlHXrXmH27J/Ttevt9O17\nNzVq6CqX4U4l7jP7968jPn641zHCQlRUNH363EW7dtcxb96veOqpZ/jud39Hp07JREdXYutW2Lq1\nYPvvfOfcEl+7FjZtKlhv3frCSvzrr79m0aJFzJ49m9mzZ5OVlcXgwYMZM2YMDzzwMrNn1z+7bdOm\nkNff32IWRYMG7WnQoD3dut3OkSPb+fjjP/Lss524+uoH6Nnz50B04MHEV1TiPuKc00i8HNSr14ZR\no95jx46lLFz4WxYuvJ+ePX9O585jqVkzeMeVO+fIzNzMli3z+OKL2UyZsoSuXbsydOhQZs2adc5N\nrwv/z+FC1anTgqFDn6Jnz5/xwQc/4bPP3uP669+gWTMdIx+OVOI+cvToTmJiqlK9ekOvo4Sl5s2v\n5NZbU9i7dw0rVjzJM88k0KhRZ+LjhxMXl8iZM525kBHt8ePHSUtLY9WqVSxevJjZsxcTFRVDy5YD\n6NHjVj755A3q1KlTbr+ehg0TGDt2EUuW/J7rruvBBx/MokePHuX2fuINlbiP5I7CI/sknwULoGrV\n4l9zLjjvERvbleHDX2LYsGfZvPlDNm2aw+rVL/D66zvp0CGe+Ph4GjduzPbtdTh+vCbO5XD6dBaZ\nmcc4enQPu3fvZuvWrezbt48OHTrQpUsXhg0bRmzsY9Su3QIzIybm3NPm69SB/v2Dk7+wqKhoEhMf\n5Ec/upwhQ4bw+uuvM3jw4OC/kXhGJe4jmko594SW8hYTU5X27UfQvv0IAJKSDlO16udkZGRw4MAB\ntmw5wuHDB4mKiiE6ujLVqtWgS5eradKkCXFxcbRp04bo6IKRe+H59dOnYf36gvXGjcunxPMNHjyC\nDh0a8v3vf58ZM2YA5fhmUqFU4j6yd+9qOnS4wesYEatWrbp069aHPn36ABAbe+7c9R135D4Xqvr1\n68c//vEPbrrpJkaPnkWdOn29jiRB4I+zHASAvXtXERur269L2SUmJvLyyy/zt7/9gCNHtnkdR4JA\nI3GfOHnyECdPHgq5k1IkOL78EqYWuptcVlb5vdfQoUMZMOBeZswYzm23/YfKlct4cLuEBI3EfWLf\nvjXExnb1zSnicmFOn4Z9+wqWwicRlYerr55EbGwP3n//dlywvhEWT6gRfGLPnlXExnb3OoaECTNj\n6NCnOHhwA2vXvuR1HLkIKnGfyJ0PV4lL8FSqVI0bb5zBggX38vnnn3sdR8pIc+I+sXfvKpKSHvI6\nhpzHCy+cu/7LX0LNmt5kCVTDhgkkJU3h1ltvZtmyT845JFL8IaCRuJlda2afmVmGmd1bzOvxZvax\nmX1jZncHP2ZkO3LkMCdOHNA1xENcTs65i1/06PETatSowRO6aasvlVrilvtN2lPAYKAjkGxm7Yts\ndgj4BfDHoCcU0tLW0LhxF6KiNEqS4DMzpk59gUceeYTNmzd7HUcuUCAj8V7AJufcdudcNjADOOcy\nes65L51zq4DT5ZAx4q1fr/nwUHD8eO6hgPnLqVNeJwqe1q1bc9999zF+/HgdreIzgcyJNwF2Flrf\nRW6xSwVZs2YFl146wusYEW/hwnOvdxJuJk2axOuvv87MmdOB0V7HkQBV+BebkydPPvs4KSmJpKSk\nio7gO6tWfcJNNz3mdQwJc9HR0fz1r3/lhz8cydix1+skIA+lpKSQkpIS0LaBlPhuoHmh9aZ5z5VJ\n4RKX0u3cuZNTp05Rp05Lr6NIBOjXrx9JSf1ZvPh/GTjwUa/jRKyiA9yHHir5yLRASjwVaGNmLYC9\nwCgg+TzbW0ApJSDLli2jR4++Z28WIP6xaRNUq+Z1igv3u989RqdOl9G16206IsoHSi1x59wZM5sI\nzCP3i9BpzrmNZnZH7svueTNrBKwEagI5ZjYJSHDOHS/P8JHgk08+oVu3Pl7HkDJ4/32vE5RNbGws\nV175G+bOvZvRoz/wOo6UIqDjxJ1zHzrn4p1zbZ1zj+Y995xz7vm8x/udc82cc3Wcc/Wcc81V4MGR\nPxIXqUi9ev2Cgwc38MUXi7yOIqXQafchLCsri3Xr1tGlS0+vo0iEiYmpwoABDzN//j3k+OnMpQik\nEg9ha9asoV27dnznO9W9jiIRqGPHH2JmLFv2d/79b84u69Z5nUwK07VTQtiyZcvo21dTKeINsygG\nDfojs2aNo1WrHxATUwWAli2hc2Tf6jWkaCQewpYsWcIVV1zhdQyJYHFxSVxySSdSU5/2OoqUQCUe\nonJycvj3v/+tk6HEcwMHPsbSpY/w9deZXkeRYqjEQ1R6ejp169aladOmXkeRCNewYQLt23+fJUse\nAeDAAfjHPwqW5cs9DhjhVOIhKiUlRaNwCRlJSZNZu/ZvHDmynRMnID29YNm1y+t0kU0lHqJU4hJK\nata8lB49fsaiRQ94HUWKUImHoPz58MTERK+jiJzVr989bNkyj3371nodRQrRIYYhKC0tjTp16mg+\nXEJKlSq1uPrqB1iw4F7GjJlb5v3Mm5d7XZl8w4eD/qiXnUbiIWju3Llcc801XscQ+Zbu3Sdw+PBW\ntmyZX+Z9HDsGBw8WLNnZQQwYgVTiIWjOnDkMGTLE6xgi3xIdXYkBAx5hwYJf45xOxw8FKvEQc+zY\nMVJTU/nud/vjHOhOWRJqOnS4gejoKnz66ZteRxE0Jx5yPvroI3r27MOf/qTrpUhoMjMGDfoj7757\nMwkJNwJVvY4U0TQSDzFz5szhmms0lSKhrUWLq2jcuDMrVuh0fK+pxEOIc47Zs2czaJBKXELfgAGP\n8J//PMbx44e9jhLRVOIhYv9++PDDFVStWp169dp7HUfCTGZm7p+x/OXMmYvfZ+7p+CN4991HLn5n\nUmaaEw8RU6fC3Ln/IDb2Jl55RffTlOB6993y2W9S0mSef/4yduyYSPPmzUv/AQk6jcRDhHOODRve\nomPHm7yOIhKwmjUvZciQX/CrX/3K6ygRSyUeInbvTiUmpgqXXHKZ11FELsiIEfeycuVK5s2b53WU\niKQSDxHp6TNJSLgJM02liL9UqVKNJ598kokTJ5KVleV1nIijEg8B2dnZrF//OpdffrPXUUTK5Lrr\nrqNDhw786U9/8jpKxFGJh4B//etf1K/flgYN4r2OIlJmTzzxBH/5y1/YvHmz11Eiiko8BEybNo2u\nXX/sdQyRi9KyZUt++9vfcuutt3ImGMcwSkBU4h7bunUrn3zyiY5KkfPaswcWLixYtmwJ7v5TUgr2\nvXhx2fdz5513EhUVxZNPPhm0bHJ+Ok7cY0888QS33347lStX18WupET5J+mUl6VLg7OfqKgoXnrp\nJXr37s2QIUNo314nrpU3jcQ9dPjwYV5//XV+8YtfeB1FJGhat27N73//e0aOHMnJkye9jhP2VOIe\nevzxxxk+fDhNmjTxOopIUE2YMIHLLruMiRMneh0l7Gk6xSP79+/nmWeeYdWqVV5HEbko27fDjBkF\n6/Hx0LWrMXXqVHr37s20adP48Y/1xX15UYl7ZMqUKYwZM4a4uDivo4hclKNHc5d89erl/rdGjRq8\n/fbbJCbFDZBzAAAGKElEQVQmEhcXx4ABA7wJGOZU4h5Yvnw577zzDunp6V5HESlX7du3Z+bMmdx0\n003Mnz+fzp07ex0p7GhOvIJ98803jB8/nj//+c/Uyx+yiISxxMREnn76aYYMGcL69eu9jhN2NBKv\nYHfddRfx8fGcPJlM4UNpdXihhKslS2Dv3pu45pocrrzyGsaM+SeXXNLT61hhQyPxCvTcc8+xcOFC\npk2bxldfGZmZnF1EwtXJk7l/xuPiRjJs2Au88spQ0tJmlP6DEhCVeAV54403mDJlCrNnz6ZWrVpe\nxxHxRHz897j55vl89NF9zJkziexsHUd+sQIqcTO71sw+M7MMM7u3hG2eNLNNZrbWzLoEN6Z/5eTk\n8Oijj/Kb3/yGDz/8kDZt2ngdSaRc5eTA6dMFS07Oua83btyFCRNWcvLkQZ599jKWLNF1yC+GuVIm\nY80sCsgABgB7gFRglHPus0LbDAEmOueGmVlv4P+cc32K2Zcr7f0qQkpKCklJSeX+Pps3b+anP/0p\nx44d46233qJp06ZnX/vzn+HYsXKPUKpt21KIi0vyOkZI0GdRoKI+i4yMD/j4418SF9eEBx54gP79\n+4fcNfUrqi/Ox8xwzhX7wQQyEu8FbHLObXfOZQMzgOFFthkOvArgnFsO1DazRheRuVylpKSU276d\ncyxfvpzbbruNPn36MHDgQJYuXXpOgYeSbdtSvI4QMvRZFKioz6Jdu+uYP38j48aNY9KkSbRt25Yp\nU6awatUqcooO4T1Snn0RDIEcndIE2FlofRe5xX6+bXbnPVeOl+zx3unTp9m3bx+bNm0iIyODZcuW\nsWjRIipXrszYsWPZtGkTdevW9TqmSEiLiYlh7Nix3HLLLaxYsYLp06fzox/9iMzMTHr06EGXLl3o\n1KkTzZo1o1mzZsTGxlKlShWvY4eMCj/EcNiwYUDuiLWs/72YnwXYsWMHc+fOLdPPZmVlceTIEb76\n6itOnjxJw4YNadu2LW3btqVXr17cc889dOjQodR/EjZqBDVqlP55lbcaNSA21usUoUGfRYGK/Cwq\nV879r5nRu3dvevfuDeT+PV29ejXr1q3jvffeY9euXezatYu9e/cSFRVFrVq1qFmzJjVq1KBSpUrE\nxMQQExNDdHT02cdRUVEl/l0839/Rwq9lZGSQmpp6QT9TkQKZE+8DTHbOXZu3/hvAOeceK7TNVGCR\nc+7veeufAYnOuf1F9uX9hLiIiA+VNCceyEg8FWhjZi2AvcAoILnINu8DPwf+nlf6R4oW+PlCiIhI\n2ZRa4s65M2Y2EZhH7heh05xzG83sjtyX3fPOudlmNtTMNgMngHHlG1tERCCA6RQREQldEX3Gppn9\nl5nlmFnEXonKzP5gZhvzTtJ628wi7nTSQE5miwRm1tTMFppZupl9amZ3ep3Ja2YWZWarzex9r7OU\nJGJL3MyaAoOA7V5n8dg8oKNzrguwCfhvj/NUqLyT2Z4CBgMdgWQzi9QbQ54G7nbOdQT6Aj+P4M8i\n3yRgg9chzidiSxz4C3CP1yG85pxb4JzLP6tiGRCaZyWVn0BOZosIzrl9zrm1eY+PAxvJPd8jIuUN\n9IYCL3qd5XwissTN7Hpgp3PuU6+zhJjbgDleh6hgxZ3MFrHFlc/M4oAuwHJvk3gqf6AX0l8chu31\nxM1sPlD41H8j9zfjt8B95E6lFH4tbJ3ns7jfOffPvG3uB7Kdc296EFFCiJnVAN4CJuWNyCOOmQ0D\n9jvn1ppZEiHcEWFb4s65QcU9b2adgDhgneWeYtUUWGVmvZxzByowYoUp6bPIZ2a3kvvPxv4VEii0\n7AaaF1pvmvdcRDKzGHIL/DXn3Cyv83ioH3C9mQ0FqgE1zexV59wtHuf6log/xNDMvgC6OecOe53F\nC2Z2LfBn4Grn3CGv81Q0M4sGPif3Kp17gRVAsnNuo6fBPGJmrwJfOufu9jpLqDCzROC/nHPXe52l\nOBE5J16EI4T/qVQB/grUAObnHUr1jNeBKpJz7gyQfzJbOjAjggu8H/AjoL+Zrcn783Ct17nk/CJ+\nJC4i4mcaiYuI+JhKXETEx1TiIiI+phIXEfExlbiIiI+pxEVEfEwlLiLiYypxEREf+39t1Ou6SgVx\n2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f096d005c10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbins = 50\n",
    "xmin, xmax = -5, 5\n",
    "\n",
    "I = quad(pdftarget, -100, +100)\n",
    "x = np.linspace(xmin, xmax, 1000)\n",
    "y = pdftarget(x, I[0])\n",
    "\n",
    "plt.hist(samples, bins=nbins, normed=True, histtype=\"stepfilled\", color=\"blue\", alpha=0.5, linewidth=0)\n",
    "plt.plot(x, y, 'k')\n",
    "plt.xlim([xmin, xmax])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make an animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_samples(samples, accept, stepsize, x, target_dist, xmin=-5, xmax=5, nbins=50, write=False, filename=\"plot_samp_mcmc.png\", trace=False):\n",
    "    nsamples = len(samples)\n",
    "    ofile = '/home/ridlo/project/stats/mcmc_sampling/'+filename \n",
    "    \n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.set_xlim(xmin, xmax)\n",
    "    \n",
    "    #x = np.linspace(xmin, xmax, 1000) \n",
    "    #target_dist = pdftarget(x, normed=normalize)\n",
    "    max_propdist = norm.pdf(0, 0, stepsize) # to draw line\n",
    "    \n",
    "    ymax = 1.1*np.amax(target_dist)\n",
    "    ax.set_ylim(0, ymax)\n",
    "    \n",
    "    ax.plot(x, target_dist, 'k') # draw target dist line\n",
    "    if nsamples > 1:\n",
    "        ax.hist(samples, normed=True, bins=nbins, histtype=\"stepfilled\", color=\"blue\", alpha=0.5, linewidth=0)\n",
    "        if trace:\n",
    "            last_state = samples[-1] # last sample \n",
    "            last_acc = accept[-1]\n",
    "            prev_state = samples[-2]\n",
    "            prev_acc = accept[-2]\n",
    "            \n",
    "            ax.plot(x, norm.pdf(x, prev_state, stepsize), 'g') # draw previous propos dist\n",
    "            ax.axvline(x=prev_state, ymin=0, ymax=(max_propdist)/(ymax), c='g')\n",
    "            \n",
    "            color = 'r'\n",
    "            if last_acc > 0: color = 'k'     \n",
    "            ax.axvline(x=last_state, ymin=0, ymax=(norm.pdf(last_state, prev_state, stepsize))/(ymax), c=color)\n",
    "\n",
    "        ratio_accept = float(len(samples[accept>0]))/float(nsamples)\n",
    "        text = r'$n_{sample} = '+'{0:d}$'.format(nsamples)+'\\n'\n",
    "        text += r'$r_{accept} = '+'{0:0.2f}$'.format(ratio_accept)\n",
    "        ax.annotate(text, xy=(0.7, 0.97), xycoords='axes fraction', ha='left', va='top') \n",
    "    \n",
    "    if write:\n",
    "        plt.savefig(ofile, bbox_inches='tight', dpi=400); plt.close()\n",
    "    else:\n",
    "        plt.show(); plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test\n",
    "mu = 0\n",
    "sigma = 2\n",
    "nsamp = 1000\n",
    "samples, accept = sampling_mcmc_mh(mu, sigma, nsamp)\n",
    "\n",
    "I = quad(pdftarget, -100, +100)\n",
    "x = np.linspace(xmin, xmax, 1000)\n",
    "y = pdftarget(x, I[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD7CAYAAACRxdTpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VFX6x/HPk4QESAIhlFAindCkVwEh9Cq4oiigNCm7\nInYWlVVAXTv7Q0VXEARFV4qoVKUIQYr0KgQTSug1hU4SMuf3R0ISMCETMpk7mTxvX/Ni7p1z731m\ngt8cztx7rhhjUEop5V48rC5AKaWU42m4K6WUG9JwV0opN6ThrpRSbkjDXSml3JCGu1JKuSEvqwu4\nSUT0nEyllMomY4xktN5lwh3A6nPux48fz/jx4y2twVWMHz+e4ODxHD+eeZv77oPOndOWIyPh22/T\nlosWheefz70anUX/XqTRzyKNK3wWIhnmOqDDMkop5ZY03JVSyg1puKcTGhpqdQkuQz+LNPpZpNHP\nIo2rfxYa7um4+g/LmfSzSKOfRRr9LNK4+mdhV7iLSBcR2S8iESIy5g7tmohIoog8lN1tlVJKOU6W\n4S4iHsBkoDNQG+grIjUyafcusCy72yqllHIse3ruTYFIY8wRY0wiMBvolUG7UcD3wNm72FYppRzC\nGMMLL7xwy7o333yThQsX8vbbb2d7XV5lT7iXA46lWz6esi6ViJQFHjTG/BeQ7GyrlFKOEhsby6RJ\nk/jtt99S1/36668A9OzZk8TERNauXWvXunXr1jn/DTiQoy5imgTkeDw9/QUBoaGhLv+FhVLKtRQr\nVoznn3+eRYsWpa5bv349DRs2BKBBgwasWrUKEbFrXatWrZz/Ju4gLCyMsLAwu9raE+4ngPLploNT\n1qXXGJgtyZdLlQC6isgNO7dNZfXVXkqprK1fv5558+YRGhqKMYa9e/fyr3/9y+qybpH+avezZ8/i\n6+sLgJ+fH6dPn8bLy8uuda7m9k7vhAkTMm1rT7hvAaqKSAXgFPAY0Dd9A2NM5ZvPRWQGsMgYs1BE\nPLPaVimVN5UrV44mTZqwYMECh+zv4sWLvPDCC0ybNi3D1/ft28eKFSsyvOR+4MCBFC1aNMPtbDYb\nnp6eACQlJeHp6Wn3urwsy3A3xiSJyNPAcpLH6KcbY8JFZETyy2bq7Ztkta3jyldKOVvLli155513\naNKkCRcvXqRw4cIO2W/hwoWpVKlSpq/XqlWLWrVq2bWv9L8AgoKCuHLlCpD8C6RUqVIAWa4rWbJk\n9t+EC7FrzN0Y8wtQ/bZ1UzJpOySrbZVSede1a9dSA33p0qV069aN9evX06JFC7777jsKFixIw4YN\n2bhxI9u3b6dnz56sXLmS3r17M2/ePLp168a+ffsYMmQIW7duZeHChXTt2pXVq1fTpEkT9u/fz8aN\nGwkKCqJr166px73Zc7+diDBgwAACAgJS16UflmnVqhVbt26la9eubN68mfbt2+Pl5cWWLVuyXJeX\nudSskEop17d3715at24NJI9NHzlyhNq1ayMi7Nu3jxIlSrBu3TomTpxIREQElStXpnTp0vj7+1Ot\nWjU8PDwoVKgQAOXLlycwMJBixYrh4+NDw4YNGT16NKNGjcLPz++W49rTc79y5QpffPEF+/fvZ9Kk\nSQwfPpx27drx888/8/333yMidOrUCWMMS5cuzXJdXiZWT7N7k4gYV6lFJZs2DZ3yV9lt3rx5hISE\ncPr0aQ4ePEjjxo1ZunQp3t7ePPLII2zdupU2bdqwdOlSSpcuTY8ePZg9ezZXr17FZrNx7do1mjVr\nRkREBPfeey82my317BWVMRHJdD53DXeVqazCvX59SH+m2JEjkO4MtLsK90uXID4+bblQIUg5gUG5\nmfj4eIYOHcpHH31EYGCg1eXkSXcKdx2WUXdt587khyOtXAm7dqUtt2wJHTs69hjKNfj4+DBr1iyr\ny3BbGu4qWy5dOsWOHV9y6dJJ7rnnPmrXfhRPzwJWl6WUuo1O+avsFhW1hs8/r8elSycoXjyE7dun\nMWPG/Vy8mOl1aUopi2jPXdnl7Nk/mDfvYR5+eA6VKrUDoFmzUaxZ8ybffNOJwYPXUqiQjpsq5Sq0\n566ylJSUwPz5fenY8YPUYAcQ8SA0dByVKnVgwYIhlt/gXCmVRsNdZWnr1s8pUiSYevUGZvh6p04f\ncOHCUXbt+trJlSmlMqPhru4oIeEya9f+mw4d3s9wTg8AT09vevSYwqpVrxIff8nJFar8zp452Pfs\n2QPAwYMHiU8513bZsmV88sknfPrpp1y7ds0ptTqThru6o507v6J8+fsJCqpzx3blyjWhcuUOrFv3\nrpMqU+qvc7VnNgd7aGgoZcuW5aeffsLHx4eYmBi+/vprRo0axdmzZ9m/f78zy3YK/UJVZcpms7Fp\n00f06vWlXe3btn2LKVPq06LFixQqFMjVq3D7hIHdu4OX/q1TDpLRXO0ZzcH+ySef0K9fv9TlOXPm\n0KxZMwDGjh2Lt7e3cwp2Iv3fTGXqwIF1eHn5cM89Le1qX7ToPdSo8Tc2bfqE0NBxJCbCjh23tkk3\nD5TKo9atW8fixYuJi4vjwoULjBw50rKbWmQ0V3tGtm7dSkBAAOHh4bz44ov88ccf+Pr6snTpUvbs\n2cOYMTm+15DL0XBXmdq48Rvq1Hk807H2jLRsOYYvv2xBixYv4e2t8wa4o5IlS+Ln50e7du1o06YN\nPj4+Dtnv5s2badq0KWD/3O32zsE+ceJERISoqCiWLVuGzWajaNGiqTNU/vzzz7fMQOkONNxVhuLj\n49m+fT7Dhu3IunE6xYtXIzj4Pv744zsaNhyaS9UpK1WvXp2tW7cyZswYChRw3NXJS5cuTQ13e+du\nv32u9ozmYJ85cyZJSUk8+eSTFCxYkN27d1O2bFnKli0LQGBgIH/88YeGu8offvnlF8qVu5eiRctn\n3fg2TZo8xa+/vkqDBk9mq9ev8gZjDAkJCanBHhUVlTp3+7hx44iJiWHDhg1cunSJzp07pz5v1apV\n6jztJUuWZOHChXTr1o24uDiMMYgIFy9epEiRInbP3Z7RXO0AR44coUKFCgCUKFEi9ZdGVFQUoaGh\n+Pj4sGrVKgBiYmKoW7durn9uzqbhrjK0cOFCGjbsfVfbVqnSiaVLR3LixGaCg5s5uDJltaNHj9Ko\nUaPU5Y8++oj//Oc/REZGUqhQId5//30+/vhjRIRRo0alPh88eHDqPO0BAQEEBgYSEBDAn3/+iaen\nJ4MGDUodP7e3557RXO1xcXH069eP9evXA9C9e3c++eQTihQpQnBwMO3aJV+It2rVKmbMmIGnpyed\n089d7SZ0yl/1FzabjbJly/LMM+tJSKhyV/tYv/4Dzp/fR69eM25Z/+qrcKcTE378UWeFzGs+++yz\n1Lnb+/Tpw6ZNm6hfvz42m43du3enPg8PD+fee+8lKSmJyMjI1HncBw4cyFdffUXNmjVp2rSpQ4d6\n3F2O53MXkS7AJNLug/reba/3BN4EbEAi8LwxZn3Ka1HAhZuvGWOaZnIMDXcXsXXrVh5//HFeemn/\nHedzv5PLl8/w6ac1eP7547d8sarhrnQed8e5U7hneRGTiHgAk4HOQG2gr4jUuK3ZSmNMPWNMA+BJ\nIP3ty21AqDGmQWbBrlzLkiVL6NGjR4724ecXRHDwffz554KsG6t85eY87hrsucueK1SbApHGmCPG\nmERgNtArfQNjzNV0i34kB/pNYudxlIv45ZdfHHLmQN26j7N79zcOqEgplV32hG454Fi65eMp624h\nIg+KSDiwCBiS7iUDrBCRLSIyLCfFqtx3+fJl9uzZQ4sWLXK8r+rVe3Hs2AYuXz7jgMqUUtnhsLNl\njDE/AT+JSCvgLeDmSGlLY8wpESlJcsiHG2MynABi/Pjxqc9DQ0MJDQ11VHnKTuvXr6dRo0apd6fP\nCW9vX6pX78nevXNo1uwZB1SnVP4WFhZGWFiYXW3tCfcTQPqTnYNT1mXIGLNORCqLSKAxJsYYcypl\n/TkR+ZHkYZ4sw11ZIywszKG/VOvWfZzVq1/XcFfKAW7v9E6YMCHTtvYMy2wBqopIBRHxBh4DFqZv\nICJV0j1vCHgbY2JEpLCI+KWs9wU6AX/Y/1aUszk63CtWbEtMzAEuXDiWdWOllMNkGe7GmCTgaWA5\nsBeYbYwJF5ERIjI8pVlvEflDRLYDnwB9UtYHAetEZAewEVhkjFnu8HehHOLKlSvs2bOH5s2bO2yf\nnp4FqF79AcLDf3DYPpVSWbNrzN0Y8wtQ/bZ1U9I9fx94P4PtDgP1c1ijcpKtW7dSp04dh4y3p1ez\nZm/Wr3+f5s2fdeh+lVKZ01MUVaqNGzemznHtSJUrd+Ts2T1cvpzxdKxKKcfTcFepNm3a5NAhmZu8\nvHyoWrUr+/f/5PB9K6UypuGugOSZ/nKr5w7JQzPh4fNzZd9Kqb/ScFcAHD9+nBs3blCxYsVc2X/V\nql04cWIzMTExubJ/pdStNNwVkDYkk1vzr3t7+1KhQhuWLfs5V/avlLqVhrsCksM9t4ZkbgoJeYC5\ncxdx8CCpj+vXc/WQSuVbGu4KuPX+lbklJKQHq1YtY+bMRGbNglmzIDo6Vw+pVL6l4a4wxrBr1y4a\nNGiQq8fx9y9DYGBVjh7NcPYJpZQDabgroqKi8Pf3p0SJErl+rJCQB4iIWJTrx1Eqv9NwV+zcuZP6\n9Z1zIXFISA8iIhahd91SKndpuCt27txJvXr1nHKs0qUbkJh4jejoCKccT6n8SsNdObXnLiKpvXel\nVO7RcFdODXfQcXelnEHDPZ+LiYkhNjaWypUrO+2YlSq14/TpnVy7plerKpVbNNzzud27d1O3bl08\nPJz3V6FAgUJUrBhKZKRerapUbtFwz+ecPSRzU7VqPYiMXOL04yqVX2i453PWhXs3Dh5cxo0bN5x+\nbKXyA7vCXUS6iMh+EYkQkTEZvN5TRHaJyA4R2SwiLe3dVlnr5rCMsxUpUo6iRSuwbdvvTj+2UvlB\nlrfZExEPYDLQHjgJbBGRBcaY/emarTTGLExpXweYC9S0c1uVC06ehMOH05YDAqB27bTlhATYuDGJ\nffv2c+FCTdavh2bNwMuuGy86RrVq3Zk9ewlBQfenrjt71nnHV8qd2fO/clMg0hhzBEBEZgO9gNSA\nNsZcTdfeD7DZu63KHceOwYoVacuVK98a7vHx8P33URQsWJL16/0BaNTIueEeEtKDhQufpFatd513\nUKXyCXuGZcoBx9ItH09ZdwsReVBEwoFFwJDsbKusce7cXkqWrJ11w1xSrlwTrlw5S1xclGU1KOWu\nHNZPM8b8BPwkIq2At4CO2d3H+PHjU5+HhoYSGhrqqPJUBs6d20fJkrUsO76IB9WqdSMiYglNm460\nrA6l8oqwsDDCwsLsamtPuJ8AyqdbDk5ZlyFjzDoRqSwigdndNn24q9x37txeKlZsZ2kN1ap1Z+fO\nGRruStnh9k7vhAkTMm1rz7DMFqCqiFQQEW/gMWBh+gYiUiXd84aAtzEmxp5tlXXOndtHqVLWDcsA\nVKnSiaNH15GQcMXSOpRyN1mGuzEmCXgaWA7sBWYbY8JFZISIDE9p1ltE/hCR7cAnQJ87bZsL70Nl\nk81m4/z5/ZQoUdPSOgoWLErZso05fHiVpXUo5W7sGnM3xvwCVL9t3ZR0z98H3rd3W2W9qKjDFC5c\nAh8ff6tLSZklcjHVqz9gdSlKuQ29QjWf+vNPa79MTa9ate5ERi7RG3go5UAa7vlUeLi1p0GmV7x4\nCF5eBTlzZpfVpSjlNjTc86n9+12n5552Aw+dSEwpR9Fwz6f+/DPc8i9T00semllsdRlKuQ0N93zI\nGMOhQ5EULx5idSmpKlRozblz+7hy5ZzVpSjlFjTc86Hz588j4kHhwsWtLiWVl5cPlSq158ABvYGH\nUo6g4Z4PRUZGUqVKNavL+IuQEL2Bh1KOouGeD7lquCffwGM5SUmJVpeiVJ6n4Z4PuWq4+/mVJjCw\nKseOrbe6FKXyPA33fCgyMpLKlataXUaGqlXrTkSEnjWjVE5puOdDrtpzBx13V8pRNNzzGWOMS4d7\nmTINuXYtlpiYg1aXolSepuGez5w5cwYfHx+KFStmdSkZSr6BR3ftvSuVQxru+UxkZCTVqrlmr/2m\nkBAdd1cqpzTc85m8EO6VK3fk+PHfiY+/ZHUpSuVZGu75TF4Idx8ff4KDm3Po0EqrS1Eqz9Jwz2fy\nQrgDOu6uVA5puOczeSXcb54SabPZrC5FqTzJrnAXkS4isl9EIkRkTAav9xORXSmPdSJSN91rUSnr\nd4jIZkcWr7LHGMPBgwfzRLgHBlbFx6cIBw7ssLoUpfKkLMNdRDyAyUBnoDbQV0Rq3NbsENDaGFMP\neAuYmu41GxBqjGlgjGnqmLLV3Th16hS+vr4UKVLE6lLsUq1aDzZt0rNmlLob9vTcmwKRxpgjxphE\nYDbQK30DY8xGY8yFlMWNQLl0L4udx1G5LK8MydwUEtKdzZt13F2pu2FP6JYDjqVbPs6t4X27oUD6\nSbkNsEJEtojIsOyXqBwlr4V7+fKtOHEiklOnTlldilJ5jpcjdyYibYHBQKt0q1saY06JSEmSQz7c\nGLMuo+3Hjx+f+jw0NJTQ0FBHlpfv5bVw9/T0pkmTrixYsIC///3vVpejlOXCwsIICwuzq6094X4C\nKJ9uOThl3S1SvkSdCnQxxsTeXG+MOZXy5zkR+ZHkYZ4sw105XmRkJH379rW6jGxp1ao38+d/ruGu\nFH/t9E6YMCHTtvYMy2wBqopIBRHxBh4DFqZvICLlgfnAE8aYg+nWFxYRv5TnvkAn4A+734lyqLzW\ncwdo3LgLmzdvJjo62upSlMpTsgx3Y0wS8DSwHNgLzDbGhIvICBEZntLsNSAQ+Oy2Ux6DgHUisoPk\nL1oXGWOWO/xdqCzZbDYOHjxI1aquOY97ZgoV8qVDhw4sWLDA6lKUylPsGnM3xvwCVL9t3ZR0z4cB\nf/my1BhzGKifwxqVA8TGniAgIAA/Pz+rS8m23r178+233zJkyBCrS1Eqz9BTFPOJM2ci81yv/aYe\nPXqwdu1aLly4kHVjpRSg4Z5vnDmT98bbbypSpAht2rRh8WK9oEkpe2m45xOnT+fdcIfkoZn58+db\nXYZSeYaGez5hT8/9xo1bH8Y4qTg79OzZk5UrV3L58mWrS1EqT3DoRUzKddkT7h9+6KRi7kJgYCDN\nmzfn559/5pFHHrG6HKVcnvbc8wGbLYlz5w7n2S9Ub+rTpw9z5syxugyl8gQN93zg4sVj+PkVp3Dh\nwlaXkiO9e/dmxYoVetaMUnbQcM8HoqMjCQrKu1+m3lSsWDHatWvHDz/8YHUpSrk8HXPPB2JiDuTZ\ncN+zB44fT1u+995+fPvtFAYPHmxdUUrlAdpzzwdiYiIpXTpvhvvFi3DkSNqjSpUebNu2jZMnT1pd\nmlIuTcM9H4iJcY9hGQBv70L06tVLv1hVKgsa7vmAu4y539S/f3++/fZbq8tQyqVpuLs5m+0GcXFR\nlCpVxepSHKZt27YcP36cP//80+pSlHJZGu5u7sKFo/j5BeHtXdDqUhzGy8uLxx57jG+++cbqUpRy\nWRrubi46OpLAQPcZkrlp8ODBzJw5k6SkJKtLUcolabi7uZgY9wz3evXqERQUxIoVK6wuRSmXpOHu\n5qKjIyle3P3CHeDJJ59k+vTpVpeRoejo5PPzbz70olrlbHaFu4h0EZH9IhIhImMyeL2fiOxKeaxL\nuVm2Xduq3OWuPXeAvn37smLFCs6dO2d1KX+xahVMm5b22L7d6opUfpNluIuIBzAZ6AzUBvqKSI3b\nmh0CWhtj6gFvAVOzsa3KRcnhnrcnDMtMQEAADzzwgH6xqlQG7Om5NwUijTFHjDGJwGygV/oGxpiN\nxpib//DcCJSzd1uVe5KSErlw4SjFilW2upRcM3ToUKZPn45xpcnnlXIB9oR7OeBYuuXjpIV3RoYC\nP9/ltsqB4uKi8Pcvi5eXj9Wl5JrWrVsTHx/P77//bnUpSrkUh04cJiJtgcFAK0fuV92dnIy3GwxX\nOcd59hNHFFc4xxXOEs9FDEnYuIHgSUGK4kNRfClJABUJoBIBVMQL5/xCERGeeuopJk+eTIsWLZxy\nTFcxfz6kv46rdWtopf/nqRT2hPsJoHy65eCUdbdI+RJ1KtDFGBObnW1vGj9+fOrz0NBQQkND7ShP\nZSY757hf5wJHWcdxfucYGzjNTgShONUpRiUKUwpfSlGU8nhSgMUygu7mc+K5SDwXOMV2wplPLIe5\nxElKUIMyNCKYZlSmAwFUzLX3OXjwYN544w1OnTpFmTJlcu04rubGDUhISFvWU/7dX1hYGGFhYXa1\ntSfctwBVRaQCcAp4DOibvoGIlAfmA08YYw5mZ9v00oe7yrmYmDufBhnDQf5kAZEs4QSbKUdT7qEl\nLfknZWiELyUz3XYxI2jMiAxfS+QaZ9jFSbYRxWpWMRZv/KlKV2rzCOVphTjwLNyAgAAee+wxpkyZ\non+HlFu7vdM7YcKETNtmGe7GmCQReRpYTvIY/XRjTLiIjEh+2UwFXgMCgc9ERIBEY0zTzLa9+7em\nsiMmJpKqVbvesi72Wixz987ly+1fsZeDVKcXzXiWSrTHG1+HHLcAhQimOcE0pykjMdg4wx4iWMRS\nnuYa0dSiDw15klLc65BjPv3003To0IFXX30Vb29vh+xTqbzMrjF3Y8wvQPXb1k1J93wYMMzebZVz\nxMQcSO257zi5hxmTJrH50nzq+naile+rdKIznhTI9ToED0pTj9LUozX/4hzh7OF/zKITgVSlMf+g\nFr3x5O5DuXbt2tSqVYt58+bRv39/B1avVN6kV6i6qcTEBC5cPE5MsUN8TQemXe+MxFVmZNIBOl+c\ni/+pHk4J9oyUpCbteJPnOEIznmE7X/DWBB828TGJXL3r/T7zzDN8/PHHelqkUmi4u62VkXPwGOLB\nMo9nqccAniOK1oylMCWsLi2VJwWoxcMMZBUAUazmIyqzng9I5Fq299ejRw9iY2NZs2aNo0tVKs/R\ncHczO07toMPXHZh87J8ERlXjKfZSjwE5GvJwlkf5kQGs5AQbmUx1dvE1Bpvd23t6evLPf/6Td955\nJxerVCpv0HB3E+eunGPEohF0/bYrD9d6mL4xL1LxQigeeewe6KW4lz7MpzffsZXPmUojogize/sn\nnniCvXv3st0NJ3PZswe2bEl7xMZmvY3KvzTc87gkWxKfbPqEWp/VolCBQoSPDOfvjf/OyWOH8vSE\nYeVpyRDWcz9j+YmB/MgArnA2y+18fHx44YUXePfdd51QpXOtWQNLlqQ9Tp+2uiLlyjTc87B95/bR\nakYrvg//njWD1jCpyySKFSoGwLFjeX+qX0GoxcM8xV58CeIz7mXVhanYzJ2HaoYPH87q1av1Nnwq\nX9Nwz4MSkhJ4c82btJnZhoH1BrJ64Gpqlax1S5tjx9xnql9v/OjEBwxgJWEXv6TjrI4ciTuSaXs/\nPz+ee+45vaBJ5Wsa7nnM3rN7afpFU34//jvbh2/n743/jofc+mO8fv06MTGnCQioYFGVuSOIuowL\nXkfHyh1p/EVjvtzxZYanPcbEwBNPPMuqVav57bdd3LhhQbFKWUzDPY8wxvDZls8I/SqUUU1HsaTf\nEu4pek+GbQ8dOkTp0hXw8MhbX6baw1O8eLnVy/w64Fc+3vQxPWf35OyVW8fiP/0UvvzSjwYNXmbY\nsNdwwXt5KJXrNNzzgHNXztFrdi+m75jOusHreLLhkyTP8pCxyMhI7rnHPYZkMlM3qC6bh23m3pL3\n0nBKQ8Kiwv7SpnHjv3P69A62bdvo/AKVspiGu4tbE7WGBlMaULNETX5/8neql8h6JofIyEiCg93z\n7kvpeXt6806Hd5jeczp95/fljTVvkGRLmxrRy6sgbdqM4403XtKrVlW+o+HuoowxTNwwkUe/f5Qv\ne33Jex3fw9vTvguRknvuIblcoevoXLUz24ZvY3XUajp904lLJu0cwfr1B3Pt2lXmzJljYYVKOZ/7\nDcq6gUvxlxiycAhRcVFsGrqJCtn8YjQiIoKHHnqE6OhcKtBChw7BpElpy0WLwuDBUNa/LCufWMmE\nNRP46HATHuF7gmmGh4cnb775Ec8805+ePXtSuHDhv+xz82bYsCFtuUoVeOABJ7wZpXKR9txdTPi5\ncJpOa0pgwUDWDl6b7WCH5HAvX949e+43bkBcXNrj4sW01zw9PHmj7Rv0kE/5jgfYwQwAmjW7nxYt\nWvDee+9luM/r12/d55UrzngnSuUuDXcXsiRiCW1mtmF0i9FMeWAKBb0KZnsfly9fJjY2llKlgnOh\nwryhhvRkML+xnndZyigSbYl88MEHfPrpp3phk8o3NNxdgDGG//v9/xi+eDgL+y5kSIMhd72vAwcO\nULVqVTw88vePtgQ1GMpm4jjMo0s6UDCwIK+//jpDhw7FZrN/MjKl8qr8nQAuICEpgRGLRzBz10x+\nf/J3mgc3z9H+IiIiqFbNvU+DtFdBitKXhQQltKTuJ/fhWaoD58/bePvt/1pdmlK5TsPdQjHXYuj8\nTWdOXT7FusHrKF+0fNYbZSEiIoKQEPccb78bggf1z79N0+uvMmZ/W6p3HsHEieM4fPiw1aUplavs\nCncR6SIi+0UkQkTGZPB6dRHZICLXReSF216LEpFdIrJDRDY7qvC8LjI6kubTmtO4TGN+evQn/H38\nHbJfDfeMNWAID/E/VgaMpvmITvTr14/ExESry1Iq12QZ7iLiAUwGOgO1gb4iUuO2ZtHAKOCDDHZh\nA0KNMQ2MMU1zWK9b2HxiM61ntmZ0i9F80OkDPD08HbZvDffMVaY9A1nN7uK/c77Oecb+a6zVJSmV\na+w5z70pEGmMOQIgIrOBXsD+mw2MMeeB8yLSI4PtBR3+SbU0cimDfhrEl72+pEdIRh9XzkRGRhIS\nEsKhQw7ftUtKSIB9+25dd6eLUUtSi2/abOQlv258uu5TWi9tTRE/x/8c8gJjIDz81nVVqoCPjzX1\nKMeyJ9zLAcfSLR8nOfDtZYAVIpIETDXGfJGNbd3KjB0zeOXXV1jYd2GOvzjNSHR0NElJSZQoUSLf\nhPuVKzDzIIUxAAAfy0lEQVR3bva2KV4wiLVD19IxqSO9v+/Np6EbgQa5Up8rs9n++tmNGqXh7i6c\ncYVqS2PMKREpSXLIhxtj1mXUMP3826GhoYSGhjqhvNxnjOHfa//Nlzu+ZM2gNXbND3M3bg7J3GlS\nMZWscIHChP0jjLb/15anNt7HCL89FC+sZxkp1xYWFkZYWJhdbe0J9xNA+tM4glPW2cUYcyrlz3Mi\n8iPJvf4sw91dJNmSeHrp02w6sYkNT26gtF/pXDuWjrdnTwHPAqx9cS21n23OlMT6jEjcSfECGvDK\ndd3e6Z0wYUKmbe0ZC98CVBWRCiLiDTwGLLxD+9Ruo4gUFhG/lOe+QCfgDzuO6Rbib8TT5/s+HIg9\nQNigsFwNdtBwt1diYvKUA9evQ3y88NEDGyh+pDqfJ9Tl1I0d2Gxprye3sbpipbIvy567MSZJRJ4G\nlpP8y2C6MSZcREYkv2ymikgQsBXwB2wi8ixQCygJ/CgiJuVY3xpjlufWm3ElVxOv8rc5f8Pf25/F\nfRfj45X7A5kRERH07t0714+T1y1alPxI48mwezfz5ZaWTK93Hwnhy4l4t3Xqq8WLJ49FK5WX2DXm\nboz5Bah+27op6Z6fATK6LdBloH5OCsyLLly/QI/velC5WGWm95yOl5PuiBQZGalXp94lDw8vBjde\nx1fr2/J10w70TvqOWt556xdlfPx1Tp+Ow2az4e/vj5+fn37/ko/plL8OFn01ms7fdKZZuWZ80u2T\nv9zfNLfYbDYN9xzy9CzAoFZr+H7TY3xf/1G6XfmUxr4jrC4rQ1eunCMqKoyoqDDOnt3DuXP7SEy8\nSGBgMUSES5cu4enpSZ06dWjfvj0PPvggDRs2tLps5UQa7g506tIpOs7qSI+QHrzT/h2n9ppOnjyJ\nv78/RYoUcdox3ZGHhyePNJ/Lij2jWVrlKaJP/0m/4v+xuiwAzp6NZN26uYSHf09s7CEqVGhNhQqh\n1KrVm5Ila9GjRxChoWl/56Kjo9m1axfLli2jd+/elC5dmrFjx9K9e3ft0ecDGu4OciTuCO2/bs+Q\nBkN49f5XnX788PBwatas6fTjuiMRoVPdDyl9sgE/FRnI+d+3M2jQIvz9HTNFRHYcOnSIuXPnMnfu\nXA4cOEmNGg/TpctH3HNPi7/cAP32vC5evDjt2rWjXbt2vP322/z000+MHj2azz77jM8//5xy5XI+\nl5FyXXrlqANEREfQemZrnmn2jCXBDhruuaFu2f4M89rKseBtBD8WzOzZs51yL9YjR47wwQcf0KRJ\nE+677z6OHj3Kf/7zH9566wTduk2mQoXWfwn2rHh6etK7d2927drF/fffT5MmTVi8eHEuvQPlCjTc\nc2j3md2EzgxlXJtxPNPsGcvq0HDPHWUK1mdsuf0Etg3kuUXP0bBRQ+bOnUtSUlLWG2fD0aNHmTRp\nEvfddx+NGjUiMjKSd999lxMnTvDZZ58RGhqKhwPmIPL29uaVV17hxx9/ZOTIv7N586cOqF65Ig33\nHNh0fBMdZ3VkUpdJObrBhiNouOeeAI9ybP7HZsq2Lkvw0GAmfTSJSpUqMWbMGHbs2HFXN/+4cOEC\nK1as4KWXXqJ27do0atSIXbt2MW7cOE6dOsXUqVNp3749Xl65M3LaokULwsLWsmnTJNaufTtXjqGs\npWPudyksKow+8/owo9cMuod0t7oc9u3bp+Gei0r6lmT1wNV0/193qv+zOgPPfsqa1fN44IE+XLwY\nQ6tWLWnUqB5Vq1YlODiYiAhfTp8uxI0b10hIuEx8/Cl+++0QBw4cYNu2bRw9epQGDRrQoUMHZsyY\nQfXqjdi6Nbln/ttvycds396x7+HSpeSbgd9kTCUGD17LjBn3U7BgMZo0+YdjD6gspeF+F27O7Djn\n4Tm0rdTW6nKIjo4mPj6esmXLWl2KWytasCjLHl/G3+b8jY9Pvs1DId8SEvI2ly6dJCRkLSdP7mPF\nihWcPHmSY8eucOnSdby8ClKggC9VqpShefNKtG3bNrW3nr5XfvYsrF176/EcHe5Xrvz1GH5+pXn8\n8eXMmHE//v5lSZ7wVbkDDfdsmrd3Hk///HSuzex4N24OyejpbbnP19uXRX0XUe/fjzKbB+nDfPz9\ny/Lgg48SFJTWbt482Ls3bblNG2hrfT8gQ8WKVeLRR3/gf//rwfPP16R4cZ3Cwh3omHs2zNw5k2d/\neZbljy93mWAHHW93Nh8vH/owj0IU4390I55LVpeUY+XKNaVt2zcZOPBvXL582epylANouNtp8ubJ\nvL76dVYPXE290vWsLodTpyAqKvmxaVM4lStruDuTJwV4kK8JJIRZdCD2eswd28fFpf28oqKSh2Fc\nTaNGw6lXrxEvvfSS3dtcvnzr+zp5MreqU9mlwzJ2eGftO0zfMZ3fBv9GxYCKVpcDwC+/wJEjyc/D\nwvbx1FMu+m9+N+aBJz34nBWMpveitqwatJwgv6AM2+7alfy4qVYt6NPHSYXaSUR4991PaNOmDsuX\nL6dTp05ZbnPwIPz4Y9pymTIwwjVnbMh3tOd+B8YYXln5Ct/s+calgv1258+HU6mS9tytIAgd+YBu\nlR6i9czWHLtwLOuNXFiRIkWZNm0aQ4cO5cKFC1aXo3JAwz0TNmPjmZ+fYfmh5awZtIay/q55JkpC\nwmWuXDlH2bKVrC4l3xKElxqPY3jD4bSe2ZoDMQesLilHOnXqROfOnXnttdesLkXlgIZ7Bm7YbjBk\nwRB2nN7BqgGrKFG4hNUlZer8+T8pXrwanp45v3pR5cyLLV7klVavEDozlGPX92a9gQt79913mTNn\nDrvSjyWpPEXH3G8TfyOe/j/051LCJZY9vgxfb1+rS7qj8+fDKVGiJkuWwNKlaeudMAWKysDwRsPx\nLeDLyIXt6cMSytLIaccOC4M1a9KW770XHnro7vZVvHhx3njjDUaOHMnatWv1NNs8SHvu6VxNvEqv\n2b2wGRsLH1vo8sEOcPbsXkqUqIkxyXezv/nQcLdO/7r9GVbmc76lK0czvl1wrrj978BdzIpwi6FD\nh3L9+nW+++47xxSonMqucBeRLiKyX0QiRGRMBq9XF5ENInJdRF7Izrau4sL1C3T5pgulfEsx95G5\nTrktniOcPbuHoKC6VpehbtO0yIM8xLfM4W8cZIXV5dwVT09PJk6cyL/+9S/i9UayeU6W4S4iHsBk\noDNQG+grIjVuaxYNjAI+uIttLRd9NZr2X7enTqk6zHxwptNui+cIZ87s1nB3UVXoyKP8yA/0Zz8L\nrC7nrrRp04aaNWsyZcqUrBsrl2JPijUFIo0xRwBEZDbJE1Dsv9nAGHMeOC8iPbK7rdWsvHtSTl2/\nHsf167EUK6Znyriq8rSiPz/zP7qTyBXq0M/ubb/99tZlq85MfOedd+jYsSODBg3SO33lIfYMy5QD\n0p+8ezxlnT1ysm2ui4qL4v4Z99O/Tn/e7fBungp2gDNn9lCyZG3ESfdpVXenLI0YwK+s4J9sY6rd\n20VG3vpITMzFIu+gbt26dO7cmYkTJ1pTgLor+TYV9p/fT+sZrXmu+XO8cv8rVpdzV3RIJu8oRW0G\nEcZa3mYDeS8kJ0yYwKeffkpsbKzVpSg72TMscwJIf7PF4JR19sjWtuPHj099HhoaSmhoqJ2HyZ4t\nJ7bQc3ZP3mn/DoPqD8qVYziDhnveEkhVBrOWWXRg3rnLPGJezzP/WqxUqRI9e/bk448/Zty4cVaX\nk2+FhYURFhZmV1t7wn0LUFVEKgCngMeAvndon/5va7a2TR/uuWXFwRX0/6E/03pOo2f1nrl+vNx0\n9uxu6tbtb3UZKhuKcg+D+I3ZcZ14bGYsz1T/Dx7iQV6YiPHVV1/lvvvu4/nnn9exd4vc3umdMGFC\npm2zHJYxxiQBTwPLgb3AbGNMuIiMEJHhACISJCLHgOeBsSJyVET8Mtv2rt9ZDn235zv6/9Cf+X3m\n5/lgt9lsnD37B6VK1bG6FJVNfgTRPzGMjUe28/fl/VmyLJ71662uKmtVq1alS5cuTJ482epSlB3s\nOufPGPMLUP22dVPSPT8D3GPvtlb4eNPHvL/+fX4d8Ct1gvJ+IJ4/H0XBggEUKlTM6lLUXShEMR5n\nGT/Qn//RnUf5AR9cvzc8duxYWrduzahRo/D397e6HHUHbv+FqjGGsb+O5dMtn7JuyDq3CHaA48d1\nvD2vK0AhHmEegVRlJqFc5ozVJWWpRo0atGvXTs97zwPcOtxv2G4wbNEwlh9azrrB61x2yt67cfz4\nbh2ScQMeeNKd/1KDB/mSFkQTaXVJWRozZgwfffQRCQkJVpei7sBtw/1q4lUenvswRy8cZdWAVZT0\nLWl1SQ4VFbWNMmWcNymVyj2C0IbXackYZtKak2y1uqQ7atCgAdWrV2f27NlWl6LuwC3D/czlM4TO\nDKWITxEW91uMv4/7jQ0eObKNsmU13N1JI4bTnf/yLV05wC9Wl3NHo0eP5sMPP8ToDHUuy+3Cfd+5\nfTSf3pxu1brx1YNf4e3pbXVJDnf69GkSEq4SEKDTDribGjzIo/zETwxiK59bXU6mOnXqhIiwYcMy\nq0tRmXCrcF91eBWhM0MZ32Y840PH55kLRLJr27ZtVKjQyG3fX35XnpYMYR0b+T+W8SI2kqwu6S9E\nhNGjR/PVVx9k3VhZIu9Mf5iFr3Z+xegVo5nz8BzaVnLvm0XfDHflOr76CtLfDOvatZztL5CqPMnv\nzKU3c+nNQ3yLN651f4FHH32U0aNf5eRJHSJ0RXm+524zNsatHseENRNYM2iN2wc7JId7xYr6P5Mr\nuXoVLl1Ke9y4kfN9FiKQx1lGQQKYSRsucTLnO3WgAgUK0L//82zYoL13V5Snw/1ywmX6zOvD8kPL\n+f3J36lZsqbVJTnF1q1bqVixsdVlKCfwxJtezKAGf2MazTnNTqtLusVDDw3l0KGVxMYetroUdZs8\nG+5RcVG0/LIlRXyKEDYwjCC/IKtLcorTp09z/fp1SpSoaHUpykkEoTVj6cgHzKIje5lrdUmpfH39\nadhwGBs3/p/Vpajb5MlwD4sKo/m05gypP4TpPafnmVviOcKmTZto0qSJfpmaD93LozzOclbwT1by\nist80dqs2TPs3v0NV69GW12KSifPhft/t/yXR79/lG8e+oZnmz+b70Juw4YNtGzZ0uoylEXK0IBh\nbOEEG/mOHlzD+vnV/f3LULPmQ2zZ8pnVpah08ky4X79xnWELhzF5y2TWD1lPh8odrC7JEhs2bKBF\nixZWl6Es5EtJHmc5xanONJpxjn1Wl8R9973Ili2fkpCQw9OElMPkiXA/FHuIFtNbcDHhIhuf3EjV\nwKpWl2SJ+Ph4duzYQdOmTa0uRVnMkwJ0YRL38yozacNe5llaT8mSNSlXrikbN35laR0qjcuf575g\n/wKGLx7Oa61fY2STkfluGCa9HTt2UK1aNZ1q1cni42H37lvXucpV9/UZRCnq8D19OMIaOjERL6z5\nDqpFi9EsWTKEpKRheKY/6V9ZwmXD/YbtBq/++ipz9s5h4WMLaRbczOqSLKdDMta4fBl++MHqKjJX\nlkYMZzsLeZIvacHDzCWQKk6vo3z5Vvj5FWfBggU89NBDTj++upVLDsscu3CMdl+1Y8/ZPWwfvl2D\nPYWGu8pMQYryCPOoxyCmcx/7mO/0GkSEjh1H8/777+uEYi7A5cJ93t55NP6iMd2qdWNJvyUUL1zc\n6pJcgjGGdevW6ZkyKlOC0IxR9GMJKxjNtJMjuZp41ak11K//INHR0axbt86px1V/ZVe4i0gXEdkv\nIhEiMiaTNh+LSKSI7BSRBunWR4nILhHZISKb73ScIQuGMHbVWBb3XczLrV7GQ1zud49TGJM8FJD+\nsWXLXgoVKkyJEhW5fBmSXOMUZ+WCytGEEWzn8o04GnzeiHUHt3P5cs7nu7GHh4cnL774Ih98oFMS\nWC3LMXcR8QAmA+2Bk8AWEVlgjNmfrk1XoIoxppqINAP+CzRPedkGhBpjsjwh10M82D5iO37efnfx\nVtzH1avw4Ye3rtu0aRXFirX/y3qlMlKQADpc/JY9fEfnWV1ozvO05J94kPtfdA4cOJBx48YRHh5O\nzZr5Y0oQV2RP17gpEGmMOWKMSQRmA71ua9ML+BrAGLMJKCoiN+cDEDuPw7Se0/J9sGfm8OFfqVy5\nvdVlqDymDn0ZzlYOsZyvCCWW3J8DplChQowcOZIPtSdiKXtCtxxwLN3y8ZR1d2pzIl0bA6wQkS0i\nMuxuC83PbLYbREWtoVKldlaXovKgopRnAL9SnQeZRlM28QkGW64e86mnnuKHH37g1KlTuXoclTln\nnArZ0hhzSkRKkhzy4caYDL9tGT9+fOrz0NBQQkNDnVCe6zt5chtFi5bH17eU1aWoPErwoAUvEkIP\nFjGUvczmAaZRktwZNilRogRPPPEEEydO1B68A4WFhREWFmZXW3vC/QRQPt1ycMq629vck1EbY8yp\nlD/PiciPJA/zZBnuKs2hQyupVEmHZFTOlaA6g1jDVj5nJq1pxnO05J94UiDD9j//DFFRacvNmkHD\nhvYda8yYMdSpU4fRo0cTFJQ/Zm3Nbbd3eidMmJBpW3uGZbYAVUWkgoh4A48BC29rsxAYACAizYE4\nY8wZESksIn4p632BTsAf9r8VBRAZuYRq1bpaXYZyE4IHTXiK4WzjGOuZSkOO8FuGbePi4MyZtMeV\nK/Yfp1y5cvTv35/333/fQZWr7Mgy3I0xScDTwHJgLzDbGBMuIiNEZHhKm6XAYRE5AEwBnkrZPAhY\nJyI7gI3AImPM8lx4H27rypWznDu3jwoV2lhdinIzRSlPP5bQhnH8wOP8QH9OXXbs3Z5efvllZsyY\nwenTpx26X5U1u85iMcb8YoypboypZox5N2XdFGPM1HRtnjbGVDXG1DPGbE9Zd9gYU98Y08AYU+fm\ntsp+ERFLqFKlE175aM565TyCUIuHGUk4RalIm9l1+XDDhyQkJThk/+XKlePxxx/X3rsF8udVQnlI\nRMRCQkIesLoM5ea88aU9/+atSr8zZ/OvVP6wDq999wNnzuR8GoGXX36ZmTNncuLE7V/Vqdyk4e7C\nEhIuc/jwaqpV62Z1KSqfOLOvGt3ilhJ67WO++PMNJsa15GjG5z/YrWzZsgwbNozXX3/dQVUqe2i4\nu7D9+xdQvnwrCuv8OsqJBKEqnRnBdhrzD36gP7PpxVn23vU+X331VRYvXszu2+dOVrlGw92F/fHH\nd9x7b1+ry1D5lOBBPZ7gaf6kPK35mnaM++NR9pzZk+19FS1alNdee42XXnpJZ4x0Eg13F3X1ajRH\nj66jRo3bZ3pQyrm8KEgLXuQZDlKjSBM6zupI77m92Xl651/axsfDgQMZP9q3H8GBA0dYtOhnC95F\n/uOyN+vI73bt+pqQkB5461w7ykV440ff8i/xYZ+nmLptKt2+7UZV34ZU4gUq0hZBiImBb77JbA8F\naNHiI559diQdO7alUKFCziw/39GeuwsyxrBt2+c0bvwPq0tR6i8KFyjMc82f4+AzB2lTpidLeZop\nNGAnM7lB/B23rVq1C/XrN+att95yUrX5l4a7CwoLW4Wnpw/33KN3XVKuq1CBQjxccThP8QcdeJc/\n+I6PqEgY47lwyzyCt3r77UlMnTqVP/7Qi9Vzk4a7C/rvfz+iSZOn8vXNwJVrMgZstrSHMclfvFal\nC4+zjCdYyRXOMYX6/I/u7Ocnkki8ZR+lS5fh7bffZsCAASQkJNyyv5v7VDknrvLNtYgYV6nFSjt3\n7qRr1248+eQhvLwKWl1OpiaIMM6BP68JE4Rx4/Tn7y4Sucpe5rGdL4jlIHUZQB36EURdnntWCAgw\nPPjgg9SoUYO6dd8jMjJt206dQG8VbB8RwRiTYS9Qv1B1MW+//TbPPPMi8fGuG+xKZaUAhanPQOoz\nkHOEs5OZzKYn3vhxY3NfhjTpy7Rp06hfvz4DBnTEx6eD1SW7HR2WcSGbNm1i/fr1DBkywupSlHKY\nktSkI+/xLIfpwVTOXj1Niy9b0P2n7nR7uxuT5z5GdExk1jtS2aLh7iJsNhujRo3inXfewc9PT39U\n7kfwoDwteT90MideOMG/2/2bgiULYvoZPvOqzeKkf3CIX0m0OWbSsvxOh2VywbJlkH6G02bNoEaN\nzNufPg2jR39GTIwnSUmPM2dO7teolJW8PLzoWKUjHat0pFnMx7w3sz8HCv/MyQZbmR/2J22OtKJd\npXa0r9SeeqXr4SHaD80uDfdccPo0HE53H+KsbgC/Z88+5s8fz5NPbuDIEf1LrPIXEeGhlt+waNFw\nYnZHMuu/ezAVtvLr4V/5YvsXRF+Npm2ltrS6pxXNg5vToEwDvD29rS7b5Wm4Wyw6OpoRI/5Gx47v\nU7x4iNXlKGUJEQ8eeGAqCxcOY+wLD7Nq1UJ61+oNwPGLx1l1eBW/H/udGTtnEBkTSb2getwXfB/N\ng5vTsExDKhWrpL3722i4W+jixYv07NmTTp0epHTpIVaXo5SlRDzo2XMap0+/TbNmzZg1axZt2rQh\nuEgwA+oNYEC9AQBcir/E1pNb+f3478zaPYsXl79I3PU46gTVoV5QveRH6XrULlkbfx9/i9+VdewK\ndxHpAkwi+QvY6caY9zJo8zHQFbgCDDLG7LR32/zm1CnYuPEEY8Y8QO3azena9R127LC6KqWsJyIM\nHDiWnj3r069fPx555BHeeOMNihQpktrG38eftpXa0rZS29R1Mddi2H1mN7tO72LTiU1M3T6V8HPh\nBBQMIKR4CNWLVyekeEjqo2JARXzc/O5mWYa7iHgAk4H2wElgi4gsMMbsT9emK1DFGFNNRJoBnwPN\n7dnWlYSFhd1yZ/HcYLPZmDZtNu+99zzNmz9P3bpj2LHD9a5EjYoKo2LFUKvLcAn6WaRx1mfRvXt3\ndu/ezYsvvki1atUYM2YMw4YNw98/4554YKFAQiuGEpquNpuxceLiCSKiI/gz+k8ioiP49fCvRERH\ncOziMQILBVK+aPnkR5Hyqc+DiwQT5BdEkG/QHX8BOCMvcsKenntTINIYcwRARGYDvYD0Ad0L+BrA\nGLNJRIqKSBBQyY5tXUZu/rCuX7/OTz/9xMSJE7l82fDYYwsIDm6eK8dyBA20NPpZpHHmZ1G8eHFm\nzpzJnj17mDBhAm+++SaPPPIIvXv3pk2bNhQseOcL/TzEg3uK3sM9Re+hfeX2t7yWZEvizJUzHL1w\nNPVxIOYAq6JWceLiCc5cOcPZK2cp5FUoNeiD/IIoVbgUJX1LUqxgMcJmh3GpzCUCCgZQrFAxihUs\nRkDBAAoXKOwSU4fYE+7l4JZZgI6THPhZtSln57Zu58aNBKKjjxATE8m5c+GsWrWGPXvW0rhxY15+\n+WWCg//Gzz/rlz9K2aNOnTp8//33nDx5klmzZqWGfKNGjahfvz716tWjUqVKlCtXjnLlytk1lbCn\nhydl/ctS1r8szTPpZBljiLsex5krZzhz+Uzqn+eunuNg7EEOxh5kyrYpxF6PJe56HLHXYom9HkuS\nLYmAggH4efvh6+2b/GcBX3y9fZP/LOCb+ppvAV8KFyiMj5cPPp4++Hj54O3pneFzH8+U5ZTnBTwL\n3PE95tYXqnf1a6t79+4AqXdqMcbc8tye13LS/tixY/zyyy93tf+EhATi4uKIjY0lISGRgIBggoJC\nKFUqhJ49+zNv3hcEBQUBsHcvlClzN5+Q8/j5ZVHjace/B1f9TLL8LPIRR3wWnp63LgcG3rpPX9+M\ntytbtixjxoxhzJgxxMbGsnnzZnbt2sXKlSs5cuQIJ06c4MSJE4gIfn5++Pn54e/vT8GCBfH09Mz0\ncbe97GsR1zCXDQEp/1WkIgA2DxuJnokkeSaR5JHEDc8bXPK4RJxHHEkeSSR5JnHD40bq8ySPJGxi\nS3543PZnRutS/jRy57mYspw4TESaA+ONMV1Sll8GTPovRkXkc2C1MWZOyvJ+oA3JwzJ33DbdPnTW\nKKWUyqacTBy2BagqIhWAU8BjwO039lwIjATmpPwyiDPGnBGR83Zse8cClVJKZV+W4W6MSRKRp4Hl\npJ3OGC4iI5JfNlONMUtFpJuIHCD5VMjBd9o2196NUkopwIXmc1dKKeU4espGJkTkRRGxiUig1bVY\nRUTeF5FwEdkpIvNFpEjWW7kPEekiIvtFJEJExlhdj1VEJFhEVonIXhHZIyLPWF2T1UTEQ0S2i8hC\nq2vJjIZ7BkQkGOgIHLG6FostB2obY+oDkcArFtfjNOkuwOsM1Ab6isgd5vZ0azeAF4wxtYH7gJH5\n+LO46Vlgn9VF3ImGe8b+DxhtdRFWM8asNMbYUhY3AsFW1uNkqRfvGWMSgZsX4OU7xpjTN6cTMcZc\nBsJJvoYlX0rp/HUDplldy51ouN9GRHoCx4wxe6yuxcUMAX62uggnyuzCvHxNRCoC9YFN1lZiqZud\nP5f+wjJfzgopIiuAoPSrSP5B/Qt4leQhmfSvua07fBZjjTGLUtqMBRKNMf+zoETlIkTED/geeDal\nB5/viEh34IwxZqeIhOLC+ZAvw90Y0zGj9SJyL1AR2CXJl60FA9tEpKkx5qwTS3SazD6Lm0RkEMn/\nBG3nlIJcxwmgfLrl4JR1+ZKIeJEc7LOMMQusrsdCLYGeItINKAT4i8jXxpgBFtf1F3oq5B2IyGGg\noTEm1uparJAyXfNEoLUxJtrqepxJRDyBP0me0fQUsBnom1+v0xCRr4HzxpgXrK7FVYhIG+BFY0xP\nq2vJiI6535nBhf/Z5QSfAH7AipTTvj6zuiBnMcYkATcvwNsLzM7Hwd4S6A+0E5EdKX8Xulhdl7oz\n7bkrpZQb0p67Ukq5IQ13pZRyQxruSinlhjTclVLKDWm4K6WUG9JwV0opN6ThrpRSbkjDXSml3ND/\nA9NFK6dQVWynAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f099b049c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_samples(samples, accept, sigma, x, y, trace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n"
     ]
    }
   ],
   "source": [
    "# make an animation\n",
    "divisor = 10\n",
    "for i in range(nsamp):\n",
    "    plot_samples(samples[0:i], \n",
    "                 accept[0:i], \n",
    "                 sigma, \n",
    "                 x, y, \n",
    "                 write=True,\n",
    "                 filename='plotsample_nsamp_{0:04d}.png'.format(i),\n",
    "                 trace=True)\n",
    "    \n",
    "    if (((i+1) % divisor) == 0):\n",
    "        print \".\","
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/zL2lg_Nfi80\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7fba40e6fe90>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "YouTubeVideo(\"zL2lg_Nfi80\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reference\n",
    "\n",
    "- SEAYAC Workshop, *MC methods in astronomy*, Tri L. Astraatmadja (MPIA Heidelberg)\n",
    "    Krabi, 4 December 2015"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}