{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PDF Sampling: The acceptance-rejection method, with different proposal distribution\n",
    "\n",
    "Acceptance-rejection method assume we have a simpler p.d.f. $q^*(x)$ which we can easily draw samples. Our proposal distribution $q^*(x)$ must cover targeted p.d.f $p^*(x)$ in our interested region.\n",
    "\n",
    "$q^*(x) > p^*(x), \\quad \\text{for all } x \\quad (\\text{or interested region } x)$\n",
    "\n",
    "$Steps :$\n",
    "- generate random variable $x$ in our interested region $[x_{min}, x_{max}]$\n",
    "- evaluate $q^*(x)$ \n",
    "- generate random variable $y$ in the range $[0, q^*(x)]$\n",
    "- if $y < p^*(x)$, then accept $x$, otherwise reject $x$.\n",
    "\n",
    "\n",
    "Lets try some proposal distribution to generate p.d.f $p^*(x)$\n",
    "\n",
    "$p^*(x) = \\exp[0.4(x-0.4)^2 - 0.08x^4]$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import quad\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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(7.8521781788662155, 1.4333237621618567e-08)\n"
     ]
    }
   ],
   "source": [
    "# integrate \n",
    "I = quad(pdftarget, -100, 100)\n",
    "print I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f309c68aed0>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUVNW59/Hv00yioAZFjM2goqABB4wiCmopRlpfFTQa\nxSHXmdeIQ5KbhdchtrlZN2pyvd7EIfEavUZN8FUBMU6I2iohjAoqNoO2IiCiOA+IDf28f+xqKJpu\n6lR3VZ0afp+1elmn6pw6T9n0r3fvs/c+5u6IiEjpqoi7ABERyS0FvYhIiVPQi4iUOAW9iEiJU9CL\niJQ4Bb2ISImLFPRmVmVmC81ssZmN28J+B5lZvZmdnOmxIiKSG5ZuHL2ZVQCLgeHAe8Bs4HR3X9jM\nfs8Aa4C73X1C1GNFRCR3orToBwNL3H2pu9cD44GRzex3KfAw8EErjhURkRyJEvSVwLKU7eXJ5zYw\ns12AUe5+B2CZHCsiIrmVrYuxtwDqfxcRKUDtI+yzAuidst0z+VyqA4HxZmbAjsCxZrYu4rEAmJkW\n3RERyZC7W5SdtvgFtAPeBPoAHYF5wN5b2P8e4ORMjw2llKbrrrsu7hJySp+vuOnzFa9kbqbN8bQt\nendfb2ZjgSmErp4/u3utmY1JnuTOpoekOzbtbx8REcmaKF03uPtTQP8mz/2phX3PS3esiIjkj2bG\n5kEikYi7hJzS5ytu+nylL+2EqXwxMy+UWkREioGZRboYqxa9iEiJU9CLiJQ4Bb2ISIlT0EtOfPop\nHHss3HQTrFkTdzUi5U1BL1n37bdw8smw884wYwbstRf87W9xVyVSviKNoxeJyh3OPx+23Rbuugva\ntYN//ANOOQX69YPvfz/uCkXKj4ZXSlbdcANMmgTPPQdbb73x+dtugyeegMcfj682kVITdXilgl6y\npqEBevYMIb/XXpu+tnZtaNGPHw+HHBJPfSKlRuPoJe/+8Q/YccfNQx6gUye45hr45S/zX5dIuVPQ\nS9Y89BCcemrLr59zDtTVwQsv5K0kEUFdN5IlDQ3Qqxc8+2zzLfpG994L990HU6fmrzaRUqWuG8mr\n6dOhW7cthzzA6afD3LmwcmV+6hIRBb1kSbpum0adOsHxx8OECbmvSUQCBb20WUMDPPxwtKCHMKb+\n4YdzW5OIbKSglzb75z9Dt83ee0fbf8QImDcPVq3KbV0iEijopc2eeAJGjYq+/1ZbwXHHqftGJF8U\n9NJms2fDkCGZHaPuG5H8iRT0ZlZlZgvNbLGZjWvm9RPNbL6ZvWJms8xsaMpr76S+ls3iJX7uMGcO\nHHhgZsdVVYXRNx98kJu6RGSjtEFvZhXArcAIYAAw2syaDqKb6u77ufsg4HzgrpTXGoCEuw9y98FZ\nqlsKRF0ddOkCPXpkdlznzmEZ44kTc1OXiGwUpUU/GFji7kvdvR4YD4xM3cHdv07Z7EII90YW8TxS\nhFrTmm900kkweXJ26xGRzUUJ4EpgWcr28uRzmzCzUWZWCzwGnJfykgPPmNlsM7uwLcVK4WlL0A8f\nDtOmQX19dmsSkU1lraXt7pPcfW9gFPDrlJeGuvsBwHHAJWY2LFvnlPi1Jeh32AH23BNmzsxuTSKy\nqSg3HlkB9E7Z7pl8rlnuPs3Mdjezbu7+sbuvTD7/oZlNJHQFTWvu2Orq6g2PE4kEiUQiQnkSl4YG\nePnltt1M5Oijw7o3w/TrXyStmpoaampqMj4u7aJmZtYOWAQMB1YCs4DR7l6bsk9fd38r+fgA4FF3\n72VmWwMV7v6lmW0DTAGud/cpzZxHi5oVmYULw3j4urrWv8fUqXDddWGJYxHJTNRFzdK26N19vZmN\nJYR0BfBnd681szHhZb8T+KGZ/Rj4FlgD/Ch5eA9gopl58lwPNBfyUpza0m3TaOhQePVV+PzzcPtB\nEck+LVMsrXbFFVBZCb/4Rdve5+ij4fLL4YQTslOXSLnQMsWSc9lo0cPGfnoRyQ0FvbTKunVhYbID\nDmj7eynoRXJLQS+tsnBh6LbZbru2v9egQfD++/Dee21/LxHZnIJeWmX+fNh//+y8V7t2cNRRatWL\n5IqCXlpl0SLo3z9773fEEfDii9l7PxHZSEEvrbJ4cXaD/vDDFfQiuaKgl1bJdot+4ED48MPQVy8i\n2aWgl4y5w5Il0K9f9t6zoiJMnprW7OIYItIWCnrJ2IoV0LVr9meyHnYYvPRSdt9TRBT00gqLFmW3\nNd9IQS+SGwp6yVi2L8Q2OvDA0CX02WfZf2+Rcqagl4xl+0Jso44dQ9hPn5799xYpZwp6ydjixbnp\nugF134jkgoJeMparFj0o6EVyQcsUS0bWrg3r23zxBXTokP33//JL2HlnWL0attoq++8vUkq0TLHk\nxFtvQZ8+uQl5gC5dYO+9Yfbs3Ly/SDlS0EtGctlt02jYMN1aUCSbFPSSkVxeiG106KEKepFsUtBL\nRvLRoh86NAyxbGjI7XlEykWkoDezKjNbaGaLzWxcM6+faGbzzewVM5tlZkOjHivFJR8t+l12Ccsr\nLFqU2/OIlIu0QW9mFcCtwAhgADDazPZqsttUd9/P3QcB5wN3ZXCsFJF8tOhhY6teRNouSot+MLDE\n3Ze6ez0wHhiZuoO7f52y2QVoiHqsFI9PPoFvv4UePXJ/rqFD1U8vki1Rgr4SWJayvTz53CbMbJSZ\n1QKPAedlcqwUh7o66NsXLO2o3bbTBVmR7GmfrTdy90nAJDMbBvwa+EGm71FdXb3hcSKRIJFIZKs8\nyYJ33oFdd83PuQYODDch+fBD6N49P+cUKXQ1NTXU1NRkfFyUoF8B9E7Z7pl8rlnuPs3Mdjezbpke\nmxr0UnjyGfTt2sGQIaGffqQ6+0SAzRvA119/faTjonTdzAb2MLM+ZtYROB2YnLqDmfVNeXwA0NHd\nP45yrBSPfAY9qJ9eJFvSBr27rwfGAlOABcB4d681szFmdlFytx+a2etm9jLwB+BHWzo2B59D8iDf\nQX/ooRp5I5INWtRMIhs4EP76V9h33/yc74svwgJnH38MnTrl55wixUSLmklWuYcWfZ8++Ttn165h\nzP7cufk7p0gpUtBLJB99FO4Atd12+T3v0KEwbVp+zylSahT0Ekm+++cbaSVLkbZT0EskcQV941II\nunwj0noKeokkrqDv2RO23jospiYiraOgl0jiCnpQ941IWynoJZI4g14XZEXaRkEvkcQd9GrRi7Se\nJkxJWu5hTPuKFfkfXgmwfj106wZvvqkFzkRSacKUZE1cY+gbtWsHhxyi5RBEWktBL2nF2W3TSN03\nIq2noJe03nkHdtst3hp0QVak9RT0klYhtOgPPhhefRXWrIm3DpFipKCXtAoh6LfZBvbZB2bOjLcO\nkWKkoJe0CiHoAQ47DF58Me4qRIqPgl7SKpSgP/xwBb1Ia2gcvWxR4xj6996DbbeNt5ZPPoHevcON\nSDp0iLcWkUKgcfSSFZ98Esaxxx3yAN/5DvTtCy+/HHclIsVFQS9btGwZ9OoVdxUbqftGJHORgt7M\nqsxsoZktNrNxzbx+hpnNT35NM7N9U157J/n8K2Y2K5vFS+4p6EWKX/t0O5hZBXArMBx4D5htZo+6\n+8KU3eqAw939MzOrAu4EhiRfawAS7v5JdkuXfCi0oD/sMLjwwrD+Tbt2cVcjUhyitOgHA0vcfam7\n1wPjgZGpO7j7DHf/LLk5A6hMedkinkcK0LJl4QJooejRA3baCV5/Pe5KRIpHlACuBJalbC9n0yBv\n6gLgyZRtB54xs9lmdmHmJUqcCq1FD6H75qWX4q5CpHik7brJhJkdCZwLDEt5eqi7rzSz7oTAr3X3\nZlctqa6u3vA4kUiQSCSyWZ60QqEG/eTJMHZs3JWI5FdNTQ01NTUZH5d2HL2ZDQGq3b0quX0l4O5+\nY5P99gUeAarc/a0W3us64At3v7mZ1zSOvgD17QtPPQV77hl3JRstWwYHHACrVkGFOgWljGVzHP1s\nYA8z62NmHYHTgclNTtabEPJnp4a8mW1tZl2Sj7cBjgHUu1okGhpg+fJwg+5C0qtXGFP/2mtxVyJS\nHNJ23bj7ejMbC0wh/GL4s7vXmtmY8LLfCVwLdANuNzMD6t19MNADmGhmnjzXA+4+JVcfRrLrgw/C\nRKnOneOuZHPDh8Nzz8F++8VdiUjh0xII0qLZs2HMmMKcifrQQ/CXv8Bjj8VdiUh8tASCtFkhXoht\nlEiEiVPr1sVdiUjhU9BLiwo56Lt3DytqzpkTdyUihU9BLy0q5KCHjf30IrJlCnppUaHNim3qqKPg\n2WfjrkKk8CnopUWF3qI//HCYNQu++SbuSkQKm4JeWvTuu4Ud9NtuCwMHwj//GXclIoVNQS/NWrcu\njKPfZZe4K9myo46CqVPjrkKksCnopVnvvRdGthT6LftGjAhLNIhIyxT00qxC759vdMghUFcH778f\ndyUihUtBL80qlqDv0CEMs5yihTVEWqSgl2YVS9ADVFXBk0+m30+kXCnopVnFFvRTpoTbC4rI5hT0\n0qxiCvqePaGyMizCJiKbU9BLs4op6AGOPVbdNyItUdBLs4ot6NVPL9IyrUcvm1m7Frp2hTVroF27\nuKuJ5ttvYaedYMmSMP5fpBxoPXpptRUr4LvfLZ6QB+jYEY48UpOnRJqjoJfNLF9eXN02jUaNgokT\n465CpPAo6GUzxdY/3+iEE8KyxV9/HXclIoUlUtCbWZWZLTSzxWY2rpnXzzCz+cmvaWa2b9RjpfAs\nWxaGLBabbt3gwAM1S1akqbRBb2YVwK3ACGAAMNrM9mqyWx1wuLvvB/wauDODY6XAFGvXDcDJJ8OE\nCXFXIVJYorToBwNL3H2pu9cD44GRqTu4+wx3/yy5OQOojHqsFJ5ibdFD6Kd//HGor4+7EpHCESXo\nK4FlKdvL2RjkzbkAaBzRnOmxUgCKuUVfWQl77AEvvBB3JSKFo30238zMjgTOBYa15vjq6uoNjxOJ\nBIlEIit1SWaK9WJso5NPDqNvjj467kpEsqumpoaampqMj0s7YcrMhgDV7l6V3L4ScHe/scl++wKP\nAFXu/lYmxyZf04SpAvDNN7DddmGyVEWRjslavBgSifCXSbF+BpEosjlhajawh5n1MbOOwOnA5CYn\n600I+bMbQz7qsVJYVqwItw8s5oDs1y+MwJk+Pe5KRApD2h9nd18PjAWmAAuA8e5ea2ZjzOyi5G7X\nAt2A283sFTObtaVjc/A5JEuK+UJsqjPPhPvvj7sKkcKgtW5kE/ffD088AX/9a9yVtM2778IBB4S/\nUDp1irsakdzQWjfSKsV+IbZR794wcKBWtBQBBb00sXx5aXTdAJx1Ftx3X9xViMRPQS+bKJUWPcAp\np8DUqfDJJ3FXIhIvBb1solQuxgJsvz0ccww89FDclYjES0EvmyjmWbHNOessjb4R0agb2WDNmtAK\nLubJUk19+234xfXii9C/f9zViGSXRt1IxlasCGvFlErIQ7jz1Pnnwx13xF2JSHxK6Eda2qqULsSm\nGjMmjL756qu4KxGJh4JeNiiloZWp+vSBYcOKfxKYSGsp6GWDUm3RA1xyCdx2G+gykJQjBb1sUEpD\nK5s6+ujQdaOFzqQcKehlg3ffDd0cpaiiAn7yk9CqFyk3CnrZ4N13wxoxpercc8ONw+vq4q5EJL8U\n9LJBKbfoIcwRuPhiuOGGuCsRyS9NmBIAPv00tOY/+wws7fSL4vXRR+HGJK+8Utp/vUh50IQpyUhj\nt00phzzADjvABRfATTfFXYlI/ijoBSj9/vlUP/95GFO/cmXclYjkh4JegPIK+p12gn/5F/XVS/lQ\n0AtQ+hdim7rqqtCqr9UdjKUMRAp6M6sys4VmttjMxjXzen8zm25m35jZz5q89o6ZzU+9abgUnqVL\ny6dFD9C9O1x9NVxxhWbLSulLG/RmVgHcCowABgCjzWyvJrt9BFwK/LaZt2gAEu4+yN0Ht7FeyZFy\n6rppdMkl4XM/9ljclYjkVpQW/WBgibsvdfd6YDwwMnUHd1/t7nOBdc0cbxHPIzEqx6Dv0AFuuQV+\n9jNYuzbuakRyJ0oAVwLLUraXJ5+LyoFnzGy2mV2YSXGSH/X1sGoV7LJL3JXk34gRMGAA/OY3cVci\nkjvt83COoe6+0sy6EwK/1t2nNbdjdXX1hseJRIJEIpGH8mTFCth559DCLUd33AGDBkFVFQwZEnc1\nIi2rqamhpqYm4+PSzow1syFAtbtXJbevBNzdb2xm3+uAL9z95hbeq8XXNTM2Pi++GEahTGv21295\nmDgR/vVfYd486No17mpEosnmzNjZwB5m1sfMOgKnA5O3dO6UIrY2sy7Jx9sAxwCvRzin5FE59s83\nddJJcOSRcNllcVcikn1pu27cfb2ZjQWmEH4x/Nnda81sTHjZ7zSzHsAcoCvQYGaXA98DugMTzcyT\n53rA3afk6sNI6yjog1tugQMOgLvvhvPOi7sakeyJ1Efv7k8B/Zs896eUx6uA5u5N9CWwf1sKlNx7\n913Yb7+4q4hfly5hqOURR4QL01VVcVckkh35uBgrBW7pUjjhhLirKAz9+8OECTByJDz9dGjhx6mh\nAebODXfGmjEDliyBdevCV2UlDB8e7p41aFDpL0gnradlioUBA2D8eNhnn7grKRwTJ8LYsfDkk7Dv\nvvk//2uvwX33wd/+BttuC4cdBoccAnvvHUZHtWsXbqAydSo89RR897tw++36HpabqBdjFfRlzj0E\nyfLlsN12cVdTWB58EC69FO6/H445Jvfnq6+HSZPg97+Ht9+Gs8+GM8+EgQO3fNz69XDXXXDtteGY\n//gP6NQp9/VK/BT0Esknn8Cuu4Ybjsjmpk2DU06BX/0KLrwwN90jn34K//M/IeB32y38chk1KvN5\nDR98ENbab2iARx5R2JcD3XhEIlm6tLxWrczUsGFhnsEf/gDHHx/+f2XLggWhe2j33eHVV+HRR8O5\nTj21dZPXdtopBHzHjuGXk5Z1kEYK+jKnoZXp9esXLogeeih8//uha+Sjj1r3Xp99BvfeG0b2/OAH\n0K3bxv74bFz47dAhXG9p3x5OOy207kUU9GWu3JYnbq2OHcOyxtOnw8KF0LcvnHsuTJmy5W6vhgZ4\n/fWwzMKoUeH/9YQJoSW/dGnoEqrMZOWoiLU++CCsXh3mBoioj77MXXEF9OoVbq8n0X34YZhY9fjj\n8PLL4TpHnz5hLH7nzqHFv3w5vPNOuE/tYYeFmbcnngjbb5+fGuvq4OCD4fnn01/QleKki7ESyYkn\nhlmgo0bFXUnxqq8PfewrV8JXX8HXX4dw79kztOB33DG+2u6+O1zknTlTF2dLkYJeIhkwIIzVjmOs\nuOSee1jHZ++9tRRzKVLQS1rusM02YS16rdhYulatCr/QZ80KI3ykdGh4paS1cmUIeIV8aevRI6zK\n+ctfxl2JxEVBX8bq6tTCKxc//WlYLmH+/LgrkTgo6MvYW2+FYYJS+rp2DcND/+3f4q5E4qCgL2Nq\n0ZeXiy6C2lp44YW4K5F8U9CXMQV9eenUCa6/PnxJeVHQlzF13ZSf00+HxYvDvXGlfCjoy5ha9OWn\nY0e45BItjVBuNI6+TH35JXTvHmZyVujXfVn5+OPwl1xtLey8c9zVSFtkdRy9mVWZ2UIzW2xm45p5\nvb+ZTTezb8zsZ5kcK/F4++2w9rlCvvx06xa6cG6/Pe5KJF/S/pibWQVwKzACGACMNrO9muz2EXAp\n8NtWHCsxeOstdduUsyuugD/9CdasibsSyYco7bnBwBJ3X+ru9cB4YGTqDu6+2t3nAusyPVbiUVen\nC7HlrH9/OOigsM6RlL4oQV8JLEvZXp58Loq2HCs5pBa9XHxxuIWhlL72cReQqrq6esPjRCJBIpGI\nrZZSV1cHxx4bdxUSpxEjwiSqBQvComdS+Gpqaqipqcn4uLSjbsxsCFDt7lXJ7SsBd/cbm9n3OuAL\nd7+5Fcdq1E0e9e8PEyfC974XdyUSp2uuCevn33xz3JVIa2Rz1M1sYA8z62NmHYHTgclbOncbjpU8\nWL8+3MZut93irkTidt554X61upF4aUsb9O6+HhgLTAEWAOPdvdbMxpjZRQBm1sPMlgE/Ba42s3fN\nrEtLx+bqw0g0K1aEOyB17hx3JRK33XcPN5159NG4K5Fc0oSpMvTss+Gm1FrcSiCMvPnf/4Wnn467\nEsmUbjwiLXrjDfXNy0YnnQRz54buPClNCvoytGCBgl422morOO00uP/+uCuRXFHQl6E33tBwOtnU\n2WeHi7LqPS1NCvoy464WvWzu4IPDaKw5c+KuRHJBQV9mPvwwhH2PHnFXIoXEDM46K7TqpfQo6MtM\n44VYS3udXsrNWWfBgw9CfX3clUi2KejLjEbcSEv69g1fU6bEXYlkm4K+zOhCrGxJ40VZKS0K+jKj\nFr1syY9+BE8+CZ9/Hnclkk0K+jKjETeyJTvsAIkETJgQdyWSTQr6MrJ6NXzzDeyyS9yVSCE76yx4\n4IG4q5BsUtCXkdpajbiR9I4/Poynf++9uCuRbFHQlxFdiJUoOncO69+MHx93JZItCvoyoguxEtWZ\nZ2rtm1KioC8jCnqJKpGAVavCvxkpfgr6MqKgl6jatYPRo3VRtlToxiNlYvXqMOvx0091MVaimTcP\nRo0KN5KvUJOwIOnGI7KJmTNh8GCFvES3336w7bbw4otxVyJtpaAvEzNmwJAhcVchxcQMzjkH7r03\n7kqkrSIFvZlVmdlCM1tsZuNa2Of3ZrbEzOaZ2aCU598xs/lm9oqZzcpW4ZKZmTPDmuMimTjzTJg0\nCb78Mu5KpC3SBr2ZVQC3AiOAAcBoM9uryT7HAn3dfU9gDHBHyssNQMLdB7n74KxVLpE1NMCsWQp6\nyVyPHjB0qJZEKHZRWvSDgSXuvtTd64HxwMgm+4wE/gLg7jOB7cys8dYWFvE8kiMLF8KOO0L37nFX\nIsVI3TfFL0oAVwLLUraXJ5/b0j4rUvZx4Bkzm21mF7a2UGk99c9LW5xwAsyfD0uXxl2JtFb7PJxj\nqLuvNLPuhMCvdfdpze1YXV294XEikSCRSOShvNI3Y4a6baT1OnWC004L69Rfc03c1ZS3mpoaampq\nMj4u7Th6MxsCVLt7VXL7SsDd/caUff4IPO/uDya3FwJHuPuqJu91HfCFu9/czHk0jj5H9tsP7roL\nDjoo7kqkWM2ZE9aqX7IkTKaSwpDNcfSzgT3MrI+ZdQROByY32Wcy8OPkiYcAn7r7KjPb2sy6JJ/f\nBjgGeD2DzyFt9MUX8OabIexFWuvAA8Na9U8/HXcl0hppg97d1wNjgSnAAmC8u9ea2Rgzuyi5zxPA\n22b2JvAn4CfJw3sA08zsFWAG8Ji7646UeTRnDuy/P3TsGHclUuwuvhjuuCP9flJ4tARCifvNb+DD\nD+HmzTrLRDLz9dfQqxfMnQu77hp3NQJaAkGSNOJGsmXrrcPNw++8M+5KJFNq0ZewdevChJdXX4XK\npgNiRVph0SI44ogw1LJTp7irEbXohZdegt12U8hL9vTvH+5S9sgjcVcimVDQl7BJk8IysyLZdPnl\n8Lvfgf4ALx4K+hLlrqCX3Dj+eFi/Hp58Mu5KJCoFfYmaNw86dNDNwCX7Kirgqqvg3/9drfpioaAv\nUY8+GlrzutGI5MIpp8DHH8Pzz8ddiUShoC9R6raRXGrXLrTqf/3ruCuRKBT0Jejtt+G99+CQQ+Ku\nRErZGWeEf2svvRR3JZKOgr4EPfoonHiiFp+S3OrQIbToL7sszNmQwqWgLzHucM89cOqpcVci5eCM\nM2D77bUGTqFT0JeYxx4LoyKOOSbuSqQcmMFtt8GvfgXvvx93NdISLYFQQtzDDUbGjYMf/jDuaqSc\njBsXrgvdd1/clZQXLYFQhqZMga++gpNOirsSKTfXXgsvvhiuD0nhUdCXCPcwgeXqq0PXjUg+dekC\nDz0EF1wACxbEXY00pUgoETU18MEH4d6eInEYPBj+8z/D/I2PP467GkmloC8Bq1bBOefATTdpSKXE\n68c/DmvhnHYarFkTdzXSSEFf5NauDX3y55yjmbBSGH77W9hxRxg+PNzdTOIXKejNrMrMFprZYjMb\n18I+vzezJWY2z8z2z+RYaR13uOgi2GUXuO66uKsRCdq3hwcegKOOCnc3q62NuyJJG/RmVgHcCowA\nBgCjzWyvJvscC/R19z2BMcAfox5bDmpqarL+nh98AKNHhwtf994b7wXYXHy+QqLPl7mKijBr9tpr\nYdiwMEjg88+zfppISv37F0WUeBgMLHH3pe5eD4wHRjbZZyTwFwB3nwlsZ2Y9Ih5b8rL5D23dujBW\ned99oXfvMKRtm22y9vatUuo/SPp8rXfOOWHJ7OXLoV8/uOEGqKvL2emaVerfvyjaR9inEliWsr2c\nEODp9qmMeKxswddfhx+M2lp4/HH4+9/DD8zf/w4HHhh3dSLp9eoV/uqcNw/++Mew2F5lZejD32ef\n8NWnD3znO1pWO1eiBH1rFOy3q74eTj45v+dctAjmzAmPUyf/ukNDQ/jv+vWhtV5fH8L9iy/gs8/C\nf3fdFfbcE0aMCGPle/XKb/0i2bD//iHob7strHg5fXq4S9VNN4UW/9dfh4u4XbvC1ltD585h4bT2\n7cOX2aZfsOkvhpZ+SSxaBHPn5v7zQahz4sT8nCsTaZdAMLMhQLW7VyW3rwTc3W9M2eePwPPu/mBy\neyFwBLBbumNT3kPrH4iIZCjKEghRWvSzgT3MrA+wEjgdGN1kn8nAJcCDyV8Mn7r7KjNbHeHYyMWK\niEjm0ga9u683s7HAFMLF2z+7e62ZjQkv+53u/oSZHWdmbwJfAedu6dicfRoREdlMwaxeKSIiuVFQ\nM2PN7FIzqzWz18zshrjryQUz+7mZNZhZt7hrySYzuyn5vZtnZo+Y2bZx19RWpTzZz8x6mtlzZrYg\n+fN2Wdw15YKZVZjZy2Y2Oe5ass3MtjOzh5I/dwvM7OCW9i2YoDezBHACsI+77wP8Lt6Kss/MegI/\nAJbGXUsOTAEGuPv+wBLg32Kup03KYLLfOuBn7j4AOAS4pMQ+X6PLgTfiLiJH/ht4wt33BvYDWuwW\nL5igBy6kwKNvAAACZUlEQVQGbnD3dQDuvjrmenLhv4BfxF1ELrj7VHdvSG7OAHrGWU8WlPRkP3d/\n393nJR9/SQiJyniryq5kw+o44K64a8m25F/Mh7n7PQDuvs7dW5x7XEhB3w843MxmmNnzZlZS04HM\n7ERgmbu/FncteXAe8GTcRbRRS5MAS46Z7QrsD8yMt5Ksa2xYleKFyN2A1WZ2T7Jr6k4z69zSzrma\nMNUsM3sG6JH6FOGbcE2ylu+4+xAzOwj4f8Du+ayvrdJ8vqsI3TaprxWVLXy+q939seQ+VwP17v7X\nGEqUDJlZF+Bh4PJky74kmNn/AVa5+7xkt3DR/byl0R44ALjE3eeY2S3AlUCzyxvmNejd/QctvWZm\n/xeYkNxvdvKC5Q7u/lHeCmyjlj6fmQ0EdgXmm5kRujXmmtlgd/8gjyW2yZa+fwBmdg7hT+Wj8lJQ\nbq0Aeqds90w+VzLMrD0h5O9z91K7CeBQ4EQzOw7oDHQ1s7+4+49jritblhN6CJJz7nkYaHHAQCF1\n3UwiGRBm1g/oUEwhvyXu/rq77+zuu7v7boRv0qBiCvl0zKyK8Gfyie6+Nu56smDDREEz60iY7Fdq\nIzfuBt5w9/+Ou5Bsc/er3L23u+9O+N49V0Ihj7uvApYlsxJgOFu46JzXFn0a9wB3m9lrwFqgZL4p\nzXBK70/JPwAdgWfCHy3McPefxFtS65X6ZD8zGwqcCbxmZq8Q/k1e5e5PxVuZZOAy4AEz6wDUkZyo\n2hxNmBIRKXGF1HUjIiI5oKAXESlxCnoRkRKnoBcRKXEKehGREqegFxEpcQp6EZESp6AXESlx/x8J\nI2t6TFu0tAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f30cc3d6fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xmin, xmax = -5, 5\n",
    "N = 100\n",
    "x = np.linspace(xmin, xmax, N)\n",
    "y = pdftarget(x, norm=I[0])\n",
    "\n",
    "plt.plot(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sampling(nsamples=1000, xmin=-5, xmax=5):\n",
    "    I = quad(pdftarget, -100, 100) # integrate\n",
    "    normalize = I[0]\n",
    "    \n",
    "    x = np.linspace(xmin, xmax, 1000)\n",
    "    p = pdftarget(x, norm=normalize)\n",
    "    yhighest = np.amax(p)\n",
    "    c = 1.01*yhighest\n",
    "    \n",
    "    # start\n",
    "    xsample = xmin + (xmax - xmin)*np.random.random(nsamples) # uniform proposal distribution\n",
    "    p = pdftarget(xsample, norm=normalize)\n",
    "    ysample = c*np.random.random(nsamples)\n",
    "    \n",
    "    accept = ysample <= p\n",
    "    \n",
    "    return xsample, ysample, accept, c, normalize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl01fWd//Hn+97cJGQhbCEQAkHEgrhhHRHF5VarRTuK\ntk6r09ppp8OxnVI9v3Zm6KI1tD3HVn+ni+O0FX9Oz3RqhzraVjujFLe0okVBlrqEVUlYQjALIQGy\nv39/3CteYiA3ZPne3Lwe53DO/S6f3FdQXvnm8/3e79fcHRERSV+hoAOIiMjgUtGLiKQ5Fb2ISJpT\n0YuIpDkVvYhImlPRi4ikuaSK3swWmtlmM9tqZktPsN/5ZtZuZh/r61gRERkc1tt19GYWArYCVwB7\ngbXATe6+uYf9ngaOAP/u7r9JdqyIiAyeZI7o5wHb3L3S3duBFcCiHvb7MvAosP8kxoqIyCBJpuin\nALsSlnfH1x1lZsXA9e7+U8D6MlZERAbXQJ2M/RGg+XcRkRSUkcQ+e4BpCcsl8XWJ/gpYYWYGTACu\nNrOOJMcCYGa66Y6ISB+5uyWz0wn/AGFgO1AKZAIbgdNPsP/PgY/1dWwsSmq76667go6QFOUcWMo5\nsJRz4MR7s9ce7/WI3t07zWwJsIrYVM9D7l5hZrfG32R59yG9je31p4+IiAyYZKZucPeVwKxu6x44\nzr5/39tYEREZOvpkbB9Eo9GgIyRFOQeWcg4s5Rx6vX5gaqiYmadKFhGR4cDMkjoZm9TUjaSfxsZG\nqqurAQiFQsycOZNQSL/giaQjFf0ItWLFU7zwQjNZWXm0te3mjjuuZdYsnUoRSUc6hBuhOjq6GDv2\nUqZOvYnMzKl0dXUFHUlEBomO6EeQJ598nvLyvwBQX99Mfv7FAScSkaGgoh9Btm7dzcGDl1FQUMqE\nCSGyswvet09XVxctLS1Hl3NycoYyoogMAhX9CJOZmc+oUWOPWdfVZTz77BrWrdtMRcVWKiuPkJWV\nA7TyjW/czMyZM4MJKyIDQkUvTJz4YV5/PXaTUbNTOfXUMwiFwlRV/YZDhw4FnE5E+ktFn+YSP5tw\nvM8pjBo1jlGjxg1VJBEZYir6NPerXz3OqlUbMTPa28NMnvyRoCOJyBBT0ae56uoDjBnzd4wde0rQ\nUUQkILqOXkQkzanoRUTSnKZu5LjC4QJ+/OPfAY8TiRhLl36KGTNmBB1LRPpIRS/HVVx8Oe5RAHbt\n+j0HDhwINpCInBQVfRrq7Oyko6MDoF/3sDEzzMJHX4vI8KSiT0M///ljrF69DbMQbW1hpkwZHXQk\nEQmQij4NvfNOE2PH3kJBwbQB+5ruYZ56ag2rV79JOGz8zd9cycSJEwfs64vI4EnqqhszW2hmm81s\nq5kt7WH7dWa2ycw2mNkrZrYgYdvOxG0DGV6GTnHx5ezbdyXbtl3Aiy9CVVVV0JFEJEm9HtGbWQi4\nH7gC2AusNbPH3X1zwm7PuPsT8f3PAh4BTo9v6wKi7t4woMllSEUiOYwffxoAhw5tCTiNiPRFMkf0\n84Bt7l7p7u3ACmBR4g7ufjhhMY9Yub/LknwfEREZBMnM0U8BdiUs7yZW/scws+uBu4FC4KMJmxx4\n2sw6geXu/uDJx5XjeeaZ1WzcuB2Aysr95Obq9IuIxAxYG7j774DfmdnFwHeBK+ObFrh7tZkVEiv8\nCndfPVDvKzFr1lRQWTmXnJwJ5OREyM8vDjqSiKSIZIp+D5B4+UZJfF2P3H21mc0ws3HuXu/u1fH1\n75jZb4n9NtBj0ZeVlR19HY1GiUajScSTd+XnFzN69JSgY4jIICkvL6e8vLzP45Ip+rXATDMrBaqB\nm4CbE3cws1PdfUf89QeBTHevN7McIOTuzWaWC1wFLDveGyUWvYiIHKv7AfCyZcet02P0WvTu3mlm\nS4BVxE6qPuTuFWZ2a2yzLwc+bmafAdqAI8An4sOLgN+amcff62F3X5X0dyUiIv2W1By9u68EZnVb\n90DC63uAe3oY9zYwt58ZRUSkH3TZo4hImlPRi4ikOV1sLSchixUryvnNb9YQCsHixddSWloadCgR\nOQ4VvfRZScmHOHJkLu3tsG/fC+zbt09FL5LCVPTSZ6FQBrm5hQBEIqMCTiMivdEcvYhImtMR/TD2\nyCP/y9q12wDYv/8QRUVZAScSkVSkoh/G3nxzFy0tHyU3t5Di4jBZWflBRxKRFKSiH+aysvLJzh4T\ndAwRSWEqeukX9xDl5evZuHEnZnDddZcxadKkoGOJSAIVvfRLcfElVFVNo6oKams3cuaZb6noRVKM\nil76JRLJobAw9tTIw4d39bK3iARBl1eKiKQ5Fb2ISJpT0YuIpDkVvYhImlPRi4ikORW9iEiaM3cP\nOgMAZuapkiWV/eEPf+TNNysBqKjYTUHB4qN3kgzanj0vYfYc4XAYM1iyZBFz5swJOpZI2jIz3N16\n2y+p6+jNbCHwI957OPj3u22/DvgO0AW0A//H3V9MZqz0zQsvvM4771xIVlYBubmZKVPyAMXFF9LZ\n+UEAqqrKqaurCziRiEASRW9mIeB+4ApgL7DWzB53980Juz3j7k/E9z8LeAQ4Pcmx0kejR09NqYJ/\nl5mRkZENxO5ZLyKpIZk5+nnANnevdPd2YAWwKHEHdz+csJhH7Mg+qbEiIjK4kin6KUDiZ9t3x9cd\nw8yuN7MK4PfA3/dlrIiIDJ4B+/3a3X8H/M7MLga+C1zZ169RVlZ29HU0GiUajQ5UPBGRYa+8vJzy\n8vI+j0um6PcA0xKWS+LreuTuq81shpmN6+vYxKIXEZFjdT8AXrZsWVLjkpm6WQvMNLNSM8sEbgKe\nSNzBzE5NeP1BINPd65MZKyIig6vXI3p37zSzJcAq3rtEssLMbo1t9uXAx83sM0AbcAT4xInGDtL3\nIiIiPUhqjt7dVwKzuq17IOH1PcA9yY6VkaGiYhsHDx4BYP788xg/fnzAiURGJl3sLINi0qTzWbNm\nFAANDTvIzX1NJ9dFAqKil0GRnV3AtGkLAOjsbAs4jcjIppuaiYikORW9iEiaU9GLiKQ5Fb2ISJpT\n0YuIpDkVvYhImlPRi4ikOV1HL0Ni7969rFu3DoDZs2eTl5cXcCKRkUPPjB0Gnn12Ndu27QZgw4a3\nKCi4lZyc4XM7gebmGurrXwHg4MFqvvCFmXz4w5cHnEpk+BvQZ8ZKsFatWk99/YVkZeWTnX3BsCp5\ngLy8IvLyrgWgsvJPQEewgURGGBX9MDF27IxhV/Aikhp0MlZEJM2p6EVE0pyKXkQkzanoRUTSnE7G\nypDbseNtOjufAeDcc89h4sTCgBOJpDcVvQypSZPmsmaNsWYNNDZWsnixc/XVVwYdSyStJTV1Y2YL\nzWyzmW01s6U9bP9bM9sU/7PazM5O2LYzvn6Dmb0ykOFl+MnKGk1p6SWUll7CmDGnBB1HZETo9Yje\nzELA/cAVwF5grZk97u6bE3Z7C7jU3RvNbCGwHJgf39YFRN29YWCji4hIMpI5op8HbHP3SndvB1YA\nixJ3cPc17t4YX1wDTEnYbEm+j4iIDIJkCngKsCtheTfHFnl3/wA8lbDswNNmttbMFvc9ooiI9MeA\nnow1sw8BnwMuTli9wN2rzayQWOFXuPvqnsaXlZUdfR2NRolGowMZT1JMKBRm9erXefPNvQB89KMX\nMnv2BwJOJZK6ysvLKS8v7/O4Xu9eaWbzgTJ3Xxhf/hrg7v79bvudDTwGLHT3Hcf5WncBTe7+gx62\n6e6Vx7F06X24fyrt7nXT1dVBY2MVAO+8U8FNN4W57rqFAacSGT6SvXtlMlM3a4GZZlZqZpnATcAT\n3d5sGrGSvyWx5M0sx8zy4q9zgauA15P/NiSdhUIZjB07QzdsExlkvU7duHunmS0BVhH7wfCQu1eY\n2a2xzb4cuBMYB/zEzAxod/d5QBHwWzPz+Hs97O6rBuubSScNDQ00NsbOb3d0tBMOBxxIRIatpObo\n3X0lMKvbugcSXi8G3nei1d3fBub2M+OI9IMf/JLdu0cRDmfQ2TmFkpL8oCOJyDClT8amqMOHOygq\nupHs7DFBRxGRYU7Xt4uIpDkVvYhImlPRi4ikORW9iEia08lYSRHGpk3bqK1tBuCyy87j1FN1d0uR\ngaCil5QwefK5VFfnUl0N9fXbKSh4U0UvMkBU9JISwuFMJk48E4D29sPAO8EGEkkjmqMXEUlzKnoR\nkTSnohcRSXOao5eUEwpl8Morm3n77dg8/ZVXns8555wRcCqR4UtFLyln0qS5NDaOY/duqKvbSknJ\nDhW9SD+o6CXlmIUYM2Y6AIcP1wF7As0jMtxpjl5EJM2p6EVE0pyKXkQkzWmOPoW8+OIr7Nq1D4BD\nh44wZkyvz/wdEZqaDrJ9+3YAiouLycnJCTiRyPBi7h50BgDMzFMlS1Buv/3/0tR0IRkZ2UQio5gw\n4XRij+AduQ4d2k9d3SrMnJaWJq67roRPfvK6oGOJpAQzw917LYmkpm7MbKGZbTazrWa2tIftf2tm\nm+J/VpvZ2cmOlWMVFZ1NcfF5FBbOGfElD5CbO5Fp0z7N1Km3kJ9/EW1tnUFHEhl2ei16MwsB9wMf\nAc4Abjaz2d12ewu41N3PAb4LLO/DWBERGUTJHNHPA7a5e6W7twMrgEWJO7j7GndvjC+uAaYkO1ZE\nRAZXMkU/BdiVsLyb94q8J/8APHWSY0VEZIAN6FU3ZvYh4HPAxSczvqys7OjraDRKNBodkFwiIumg\nvLyc8vLyPo9Lpuj3ANMSlkvo4TPp8ROwy4GF7t7Ql7HvSix6ke7C4Sz+/OfNbNnybwAsXHg+F188\nL+BUIkOn+wHwsmXLkhqXzNTNWmCmmZWaWSZwE/BE4g5mNg14DLjF3Xf0ZaxIsiZMmE1u7j/Q2vo3\n7N07l+3bd/U+SER6P6J3904zWwKsIvaD4SF3rzCzW2ObfTlwJzAO+InFrglsd/d5xxs7aN+NpDUz\nIze3EIDm5hpgX7CBRIaJpObo3X0lMKvbugcSXi8GFic7VkREho7udSMikuZU9CIiaU5FLyKS5lT0\nIiJpTkUvIpLmdD96GbZaW1upr68HID8/n0gkEnAikdSkog9YdXX10bLq7OwIOM3wkZMzgT//uY5X\nXvklHR1tzJ9fyLx5ZwJQWlrKhAkTAk4okjr04JGAfeUr99LQUEIoFMY9h5KSqwmFwkHHGlba2pqp\nri7HrIuWlkYuvTSDL3zh5qBjiQy6ZB88oiP6gLW3O8XFi4hE9Hi8k5WZmUdp6V8DUFu7ha6u9QEn\nEkktOhkrIpLmVPQiImlORS8ikuZU9CIiaU4nYyWthEIZvPlmJXff/RAACxacyaWXXhBwKpFgqegl\nrYwdO4Ompk+zf38XBw/upr39ZdxbAJg+fTqlpaUBJxQZeip6SStmxujRJQDk5U3i7bfb2bGjk5aW\nA8ybt5OvfvXvAk4oMvRU9JK2wuFMSksvA6Ch4S3ghWADiQREJ2NFRNKcil5EJM0lVfRmttDMNpvZ\nVjNb2sP2WWb2kpm1mNlXum3baWabzGyDmb0yUMFFRCQ5vc7Rm1kIuB+4AtgLrDWzx919c8JudcCX\nget7+BJdQNTdGwYgr8hJiURy2LRpF1/84t0ARKNn8clP/nXAqUSGRjInY+cB29y9EsDMVgCLgKNF\n7+61QK2Z9fQvx9AUkQQsL28So0b9C+5dNDbuorJyddCRRIZMMgU8BdiVsLw7vi5ZDjxtZmvNbHFf\nwokMpHA4k4yMbMLhzKCjiAypobi8coG7V5tZIbHCr3D3Hg+nysrKjr6ORqNEo9EhiDf0tmzZQl1d\nHQAdHe0BpxGR4aK8vJzy8vI+j+v1wSNmNh8oc/eF8eWvAe7u3+9h37uAJnf/wXG+1nG3j6QHj/zj\nP97NoUPnEAqFCYdzmTJlAWa9PjtABsiBA5WMGfME119/CQCTJ0+mqKgo4FQifTeQDx5ZC8w0s1Kg\nGrgJONHje46+qZnlACF3bzazXOAqYFkS75nW3GHq1CvIyMgKOsqIlJdXRGXlKfzrv+6kre0QZ575\nKkuXfj7oWCKDpteid/dOM1sCrCI2p/+Qu1eY2a2xzb7czIqAdUA+0GVmtwNzgELgt2bm8fd62N1X\nDdY3I5KMjIzso0+kOnhwN52dKwNOJDK4kpqjd/eVwKxu6x5IeF0DTO1haDMwtz8BRQZbS8sRtmzZ\nAkBhYSHjxo0LONHQ2bt3L7t37wYgHA4zd+5cwmE9szjd6F43MqKNGjWeXbsmc++962lvb2H27A6+\n/vX0uDjswQd/zcsvbwXArJPS0klEIpmMGpXB5z//cXJzc3n44T+wYcMosrNH09DwGmPG/J5wOEw4\n3Mn111/AhAkTCIVCnHPOOWRkqC6GK/2XkxEtEhlFaemNADQ1VdPc/N/s3LnzffuNHTuWgoKCIU7X\nP7t311NQ8Hlyc4toa2umrq4RgNra59i27X5CoTAHDrRQXLyYvLwipk69GvcuABoatvLII9uBGlpb\nt3PXXXnMmjXrBO8mqUxFLxKXnT2GPXvG853vPH/M+s7Odk49tZM77/xiQMlOnlmYUChMdnYB2dmx\nH1R5eZ+io6Ml/jpMJJIT39cwi03bjB9/OuPHnw7Arl3/xbtXxHV1ddHe3n50/8xMfSZhOFDRi8RF\nIqOYPv1T71vf0tLIjh33U1b2MwDmzJnKJz7x0SHJVFNTw4EDBwDIyclh6tTYqbDt27ezf/9+ALKz\ns5k7dy5mxp49e6ipqQHgyJHDPX7NcDhCOBzpU47NmzdTX1/Piy+uY9OmeiKRTEKhDr75zU8zY8aM\nk/32ZIio6EV6kZ1dQFHRF2lubqOjo4XHH/8vnn12EwBnnTWVJUtuGbT3/t73/oPGxmJCoTCwg/vu\n+wo5OTn89KdPsH//aWRkZFFbu4px457AzDh82MnKmkMkko3ZmRQX9//Ecn7+fB5/fBvQjPvplJYu\nIBLJoarqMZqbm/v99WXwqehFkjBq1HuFmZ//T7h30dbWxPbtP6e2thaASCRydCojIyODSKRvR809\naW3tYvLkjxOJjKKq6l5eeuklMjMzaW1tZcqUS8nOLqC09Kp+v8+JjBlzCmPGnDKo7yGDS0Uv0kfv\nTnuYGfv2jWPp0hWA09raSGZmBHentDSXsrIvn9TXf/XVDdTWxm6R8e58OEBW1kJ++ct3gA4yMqKM\nH5/f329FRggVvchJCoczmTHj/Z+o7ehoYf/+H/U45vDhw7S2tgLwwx/+B1VVsfn3kpICvv3t2wmF\nQixf/hStrRcRCmWQmXktGRnZABQVnTVI38nJcc/gN7/5IytXvkpGhvHpT1/DxIkTg44lPVDRiwyi\nzs5O7rvvF+zb1wQ4NTUNhEJjAOjoGMMpp9wOQGXld6ipqSEUCuHexdSpF6X8XTaLi6/iwIH9HDgA\ntbUvcPnle1X0KUpFLzIIurq6aGhooL29nU2b9lBU9I8ATJmSffRyxmN9gDvu+C0AHR2lRy9zTGWR\nyCjGjCkFoKlpQ8Bp5ERU9EOkra3t6HzrSLlL50gVCmVw+HAR//RPv4ivOeWYk7k9mT79psEPJiOW\nin6ILFv2b+zd246Z0dY2lokTU/+ITU5OKJTR49y9SFBU9ENk//5DlJQs7fMHVURE+kvPchURSXM6\noheRfnOHjRtfp7r6HQAuvng++fm6zj9VqOhFpN/Gj1/An/4UuyVyU9ObTJw4jvPOOy/gVPIuFb2I\n9FtubiG5uYUAVFXV09nZSUdHB4DuY58C9F9ARAZUJDKOn/50FbGnj3byL//yCU4//fSgY41oKnoR\nGVCTJy8AFgBQVfUkTU1NwQaS5K66MbOFZrbZzLaa2dIets8ys5fMrMXMvtKXsSIiMrh6LXozCwH3\nAx8BzgBuNrPZ3XarA74M3HsSY0VEZBAlc0Q/D9jm7pXu3g6sABYl7uDute7+KtDR17EiIjK4kpmj\nnwLsSljeTazAk9GfsSIyzLmHePbZdaxbtwMzuOGGDzFp0qSgY404KXUytqys7OjraDRKNBoNLIuI\n9F9x8aXs2TOdPXugtnYD06dv4IwzzgBgypQphMO651NflJeXU15e3udxyRT9HmBawnJJfF0y+jQ2\nsehFZPiLRHKYMCF2Ws4MHn30RR59dC8tLQe47baLueCCCwJOOLx0PwBetmxZUuOSKfq1wEwzKwWq\ngZuAm0+wv/VjrIikqfHjZzN+fKz0d+58+pjHJMrg6rXo3b3TzJYQ+/RDCHjI3SvM7NbYZl9uZkXA\nOiAf6DKz24E57t7c09hB+25EROR9kpqjd/eVwKxu6x5IeF0DTE12rIiIDJ2UOhmbbr7znX+jqqoB\ngObmLGIfKxARGVoq+kH01lu1lJR8HTDMQoRCusJARIaein6QhUIZOpIX6UF1dTWbNm0CYNasWWRn\nZwecKH2p6EVkyI0bdzorV65l5cq3OHiwmttuO8gll1wSdKy0paIXkSE3enQJo0eXALBz53O4e8CJ\n0pvmFERE0pyKXkQkzWnqRkQC5Q6vv76V+vpDAFx00V8xcWJhwKnSi4peRAJVXHw+69fnsH491Ndv\nZ/To17n88g8FHSutqOhFJFBZWfmUlMwHoKOjFegKNlAa0hy9iEia0xH9AGtoaODw4cNBxxAROUpF\nP4C6urq4886f0NIyMb7mdI69a7OIyNBT0Q+w5uYOpk9fHHQMkWGrra2N5uZmAHJzczHTwVJ/qehF\nJGXk5hbyyCNP8sgjr9HR0cbixVEuueSioGMNeyp6EUkZhYVzgDkAVFWt1vmuAaKrbkRE0pyO6EUk\nZdXV1bFlyxYASktLdSvjk6SiF5GUNHbsDJ59dhfPPbeeQ4fqueWWD3DNNVcGHWtYSmrqxswWmtlm\nM9tqZkuPs899ZrbNzDaa2bkJ63ea2SYz22BmrwxUcBFJb/n5xUydejMlJTeTk3MuXV36xOzJ6vWI\n3mKPR7ofuALYC6w1s8fdfXPCPlcDp7r7aWZ2AfBTYH58cxcQdfeGAU8vIiK9SuaIfh6wzd0r3b0d\nWAEs6rbPIuAXAO7+MlBgZkXxbZbk+4iIyCBIZo5+CrArYXk3sfI/0T574utqAAeeNrNOYLm7P3jy\ncVPTG2+8wcGDB/WUHBFJSUNxMnaBu1ebWSGxwq9w99U97VhWVnb0dTQaJRqNDkG8/uno6ODeex8l\nNmMFmZkLA04kkn4yMrJZtepPvPTSNgA+/vFLOO+8cwJONfTKy8spLy/v8zjr7SjUzOYDZe6+ML78\nNcDd/fsJ+/wMeN7dfx1f3gxc5u413b7WXUCTu/+gh/fx4XhE3NHRweLF36O09I6go4ikLXfnyJE6\nAGpq/sKiRS187GPXBJwqeGaGu/d6j4hk5s7XAjPNrNTMMoGbgCe67fME8Jn4G88HDrh7jZnlmFle\nfH0ucBXweh++DxERzIycnAnk5EwgEskNOs6w0+vUjbt3mtkSYBWxHwwPuXuFmd0a2+zL3f1JM7vG\nzLYDh4DPxYcXAb81M4+/18PuvmpwvhUREelJUnP07r4SmNVt3QPdlpf0MO5tYG5/AoqISP/oskcR\nkTSnWyCcpJqaGg4fPkxnZ2fQUURGnKamg7z99tsATJ48WffA6UWvV90MleF01U1bWxu33XYP7e0l\n8TXjmTr12kAziYwUTU17aWhYhRm0tDRx440zueGGq4OOFYhkr7rREf1JcHdaW8OUln426CgiI05+\nfjH5+Z8FYO/edbS3VwcbaBjQHL2ISJrTEb2IDGs7dlTx2GP/C8AHP3gmp5xSGnCi1KOiF5Fha+LE\ns3jrLeett+DgwV20tm5U0fdARS8iw1ZGRhZTppwPQCiUAVQFGyhFqeiT1NbWxre+dR81NYcAaG8v\nDDiRiEhyVPRJ6ujooKamk2nT7oyv6fWKJhGRlKCi7wMzI/bALRFJNWZGRcXbPPTQfwNw4YVnM2fO\nrF5GjQwqehFJCxMnnkltbYS6OqexsYqMjDdU9HEqehFJC6FQBhMnngGAexewPdhAKURFfwJtbW3c\nffdyamsP4+60teUFHUlEkhCJ5PDHP77BunWxsr/22gu46qrLAk4VHBX9CbS2trJz5xEmT/4SAGPH\nZgacSESSMW7cTPLzv4q7U1e3haqq7cc809lsZF1MoaLvhZmRmakn2ogMN5FIDgA5OYU8//z/sHr1\nt3F3rr9+PjfcMLKe7ayi74G7H/0jIsNbQcFUCgq+BcA771RQU7Mp4ERDT0XfzZEjR/jmN39MfX0r\nAF1dJb2MEJHhwizE9u27+dnPVgBw/vlzOO+8swNONfiSKnozWwj8iPeeGfv9Hva5D7ia2DNjP+vu\nG5MdGzR3Z/369bS0tNDS0sKBA1lMn/61oGOJyAAbP/406uqu5bXXumhq2ktr60bmzDkNgMzMTMLh\ncMAJB0evn/6x2CeE7gc+ApwB3Gxms7vtczVwqrufBtwK/CzZsamgvr6eH//4aR588BD/+Z+dZGdf\n0eN+O3eWD22wk6ScA0s5B1aQOc1CTJgwi8LC0ykqOot162r40pfu44tf/CEPP/z4MfuWl5cHE3IQ\nJHNEPw/Y5u6VAGa2AlgEbE7YZxHwCwB3f9nMCsysCDglibGB6Ozs5PnnX6ClpZWWlhbC4RxKS686\n4ZidO8uZPj06NAH7QTkHlnIOrFTJmZs7kVNP/WcAmpv3sWrVg7zwwpsAXHHFOWze/CrRaDTAhAMn\nmaKfAuxKWN5NrPx722dKkmMHlbuzc+dO2tvbaW9vp6LiLdyNlpYjrFxZSXb2RUA+o0efNZSxRCSF\n5OVNYsaMrwPQ1FTNM8/8gq1bX+HIkR9SWlrA7NmnADB16lROO20mAI2NjXR0dACQm5ub0s+tHayT\nsUNykerBgwf53vd+QFfX+7e1tsK7z+1ubIS2ttjrI0cg/t+GzMwcsrN3AtDcvJPm5jUnfL/GxtfY\ntetXA5R+8CjnwFLOgTUcch46VEBdXRc7dnTxxhtVPPlk7PbHxcWQE7tq85guGT0asrJir9vb3+ue\nnBzIjH/8JrGnzj77bG688WND8J3E9PpwcDObD5S5+8L48tcATzypamY/A55391/HlzcDlxGbujnh\n2ISvoWtIKGt5AAAD50lEQVQZRUT6aKAeDr4WmGlmpUA1cBNwc7d9ngC+BPw6/oPhgLvXmFltEmOT\nDisiIn3Xa9G7e6eZLQFW8d4lkhVmdmtssy939yfN7Boz207s8srPnWjsoH03IiLyPr1O3YiIyPCW\nck/RMLOvmlmXmY0LOktPzOzbZrbJzDaY2UozmxR0pp6Y2T1mVmFmG83sMTMbHXSmnpjZjWb2upl1\nmtkHg86TyMwWmtlmM9tqZkuDznM8ZvaQmdWY2V+CznI8ZlZiZs+Z2Rtm9pqZ3RZ0pp6YWZaZvRz/\n9/2amd0VdKYTMbOQma03sydOtF9KFb2ZlQBXApVBZzmBe9z9HHc/F/hfIFX/R1gFnOHuc4FtwNcD\nznM8rwE3AH8MOkii4fJhv7ifE8uZyjqAr7j7GcCFwJdS8e/T3VuBD8X/fc8FrjazIb0kvI9uB97s\nbaeUKnrgh8A/Bx3iRNy9OWExF+jh4s7gufszHnv6AsAaICVv2uPuW9x9G6n3EN6jHxR093bg3Q/7\npRx3Xw00BJ3jRNx937u3RYn/G6og9jmblOPuh+Mvs4idx0zJ+e34gfE1wP/rbd+UKXozuw7Y5e6v\nBZ2lN2b2XTOrAv4W+FbQeZLw98BTQYcYZo73IUDpJzObTuxo+eVgk/QsPh2yAdgHPO3ua4POdBzv\nHhj3+oNoSO9eaWZPA0WJq4iFvAP4BrFpm8RtgThBzm+6++/d/Q7gjvi87ZeBsqFP2XvO+D7fBNrd\nPbBPqCSTU0YGM8sDHgVu7/bbccqI/yZ8bvy81u/MbI679zo9MpTM7KNAjbtvNLMovfTlkBa9u1/Z\n03ozOxOYDmyy2KNfSoBXzWyeu+8fwojA8XP24FfAkwRU9L3lNLPPEvvV7vIhCXQcffj7TCV7gGkJ\nyyXxdXKSzCyDWMn/p7s/3tv+QXP3g2b2PLCQJObBh9gC4DozuwYYBeSb2S/c/TM97ZwSUzfu/rq7\nT3L3Ge5+CrFfk88NouR7Y2YzExavJzbXmHLit4f+Z+C6+Amm4SCV5umPflDQzDKJfdjvhFc2BMxI\nrb+/nvw78Ka7/zjoIMdjZhPMrCD+ehSxWYbAb8LYnbt/w92nufsMYv9vPne8kocUKfoeOKn7P+33\nzOwvZrYR+DCxs96p6F+BPODp+OVXPwk6UE/M7Hoz2wXMB/7HzFLiXIK7dwLvftjvDWBFqn7Yz8x+\nBbwEfMDMqszsc0Fn6s7MFgCfAi6PX7q4Pn4wkmomA8/H/32/DPzB3Z8MOFO/6QNTIiJpLlWP6EVE\nZICo6EVE0pyKXkQkzanoRUTSnIpeRCTNqehFRNKcil5EJM2p6EVE0tz/BzlBhUJ6uyQQAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f30a95c59d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xsample, ysample, accept, c, normalize = sampling(nsamples=1000000)\n",
    "xaccept = xsample[accept]\n",
    "plt.hist(xaccept, bins=100, normed=1, histtype='stepfilled', color='blue', alpha=0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make an animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_samples(xsample, ysample, accept, c, normalize, xmin=-5, xmax=5, nbins=50, write=False, filename='plot_rejection_sampling.png', trace=False):\n",
    "    nsamples = len(xsample)\n",
    "    ofile = '/home/ridlo/project/stats/rejection_uniform/'+filename\n",
    "    \n",
    "    xaccept = xsample[accept]\n",
    "    n_accept = len(xaccept)\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.set_xlim(xmin, xmax)\n",
    "    ymax = 1.2*c\n",
    "    ax.set_ylim(0, ymax)\n",
    "    \n",
    "    x = np.linspace(xmin, xmax, 100) # draw target dist line\n",
    "    target_dist = pdftarget(x, norm=normalize)\n",
    "    propos_dist = np.empty_like(x)  # draw uniform dist line\n",
    "    propos_dist.fill(c)\n",
    "    \n",
    "    ax.plot(x, target_dist, 'k')\n",
    "    ax.plot(x, propos_dist, 'g')\n",
    "    \n",
    "    if nsamples > 0:\n",
    "        if n_accept > 0:\n",
    "            ax.hist(xaccept, bins=nbins, normed=1, histtype='stepfilled', color='b', linewidth=0, alpha=0.5)\n",
    "    \n",
    "            if trace:\n",
    "                sample_x = xsample[-1] # last sample \n",
    "                sample_y = ysample[-1]\n",
    "                sample_accept = accept[-1]\n",
    "                ax.axvline(x=sample_x, ymin=0, ymax=(sample_y)/(ymax), c='k')\n",
    "                ax.axhline(y=sample_y, xmin=xmin, xmax=(sample_x-xmin)/(xmax-xmin), c='k')\n",
    "                if sample_accept:\n",
    "                    ax.plot(sample_x, sample_y, 'k.')\n",
    "                else:\n",
    "                    ax.plot(sample_x, sample_y, 'r.')\n",
    "\n",
    "            rate_accept = float(n_accept)/float(nsamples)\n",
    "            text = r'$n_{sample} = '+'{0:d}$'.format(nsamples)+'\\n'\n",
    "            text += r'$r_{accept} = '+'{0:0.2f}$'.format(rate_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=300); plt.close()\n",
    "    else:\n",
    "        plt.show(); plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAD7CAYAAACc26SuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8FHUe//HXNxVCT2hHCy0JBEIAaVIDFspBACVKOymH\n9aeex52i553AwenvVKR49vZTUUB68ADFEoQgJsABCWkI0iEQekgCyeb7+yM9pCxhs7Oz+3k+Hnnc\nzuxk9p3Fe2fy3ZnvKK01QgghzMnN6ABCCCGqTkpcCCFMTEpcCCFMTEpcCCFMTEpcCCFMTEpcCCFM\nzMOeL6aUkvMZhRCiCrTWqqz1dj8S11ob/jV79mzDMzjKl7wX8l7Ie+H470VFZDhFCCFMTEpcCCFM\nzCVLPCwszOgIDkPeiyLyXhSR96KIo78XqrLxFpu+mFLanq/nCj75BC5eLFp+4gmoUcO67z17FpYu\nLf/5jh1h+PDbyyeEuH1KKXQ5H2za9ewUYXvp6XDlStHyrfyOzM0t+b2lZWZWPZcQwj5ccjhFCCGc\nhZS4EMKUtNbMnDmzxLp58+YRGRnJyy+/fMvrzEpKXAhhOhcvXmTRokX89NNPheu+//57AMLDw8nO\nzmbbtm1Wrdu+fbv9fwAbkhIXQphOgwYN+POf/0zdunUL10VHR9OtWzcAunXrxg8//GD1OjOTDzaF\nEGWKjo5m5cqVhIWFobXmwIED/P3vfzc6VgnFz3Y7e/YstWrVAqB27dqcOXMGDw8Pq9aZmZS4EKJC\nzZs3p2fPnqxfv94m+7ty5QozZ87kww8/LPP5hIQEtmzZglI3n1E3ZcoU6tWrV+b35ebm4u7uDoDF\nYsHd3d3qdWYmJS6EKFO/fv145ZVX6NmzJ1euXMHHx8cm+/Xx8aFNmzblPh8cHExwcLBV+ype9E2a\nNOHatWtA3i+Kxo0bA1S6rlGjRrf+QzgQKXEhRJkyMzMLi3vjxo2MGDGC6Oho+vbty7Jly6hRowbd\nu3dn586d7Nmzh/DwcL777jvuv/9+Vq5cyYgRI0hISGD69Ons2rWLyMhIhg8fzo8//kjPnj1JSkpi\n586dNGnShOHFriorOBIvTSnFQw89RP369QvXFR9O6d+/P7t27WL48OHExMRw11134eHhQWxsbKXr\nzExKXAhRpgMHDjBw4EAgb+z46NGjdOrUCaUUCQkJNGzYkO3bt7NgwQJSUlJo27YtTZs2pU6dOgQE\nBODm5kbNmjUBaNWqFb6+vjRo0ABvb2+6d+/Os88+y1NPPUXt2rVLvK41R+LXrl3jgw8+ICkpiUWL\nFvHII48wZMgQNm3axKpVq1BKce+996K1ZuPGjZWuMzO57N7k3nwTzp8vWp41C/L/f1OpM2fg3XfL\nf75LF7jvvtvLJ5zPypUrCQwM5MyZMxw6dIgePXqwceNGvLy8iIiIYNeuXQwaNIiNGzfStGlTRo4c\nyfLly8nIyCA3N5fMzEx69+5NSkoKnTt3Jjc3l+7duxv9Yzm0ii67lxI3OSlx4eiuX7/OjBkzWLx4\nMb6+vkbHMSWZO0UIYRhvb28+//xzo2M4LSlxUWWnTuXNhFigZUvw8zMujxCuSEpcVFlCAhS/Yjk8\nXEpcCHuTy+6FEMLEpMSFEMLEpMSFEMLEZExcCGEK8+bNIzQ0lPj4eP72t7/d9Hx2djbvv/8+WVlZ\nXLp0iXnz5pGbm8uyZcuoWbMmqampPP744wYkr15yJC6EcHjWzAG+atUqJk6cyF/+8heSkpKIiYlh\n8+bNhISEcN9999GkSRP27t1r7+jVzqoSV0oNU0olKaVSlFKzKtiup1IqWykll4gIIWzGmjnAk5OT\nWbFiBQBt27blxIkT1KlTh5deeolr165x6tSpCifeMqtKS1wp5Qb8BxgKdAImKKU6lLPd/wW+sXVI\nIYT9bd++neeff57HHnuMCRMmGHoHnLLmCi/thRdeYMqUKQDs37+f3r17M2DAAHx9fenUqRO1a9cu\ndxpbM7NmTLwXcFBrfRRAKbUcGA0kldruKWAV0NOmCYUQhmjUqBG1a9dmyJAhDBo0CG9vb5vsNyYm\nhl69egHWzx1uzRzgBfm2b9/OkCFDaN68OWfOnKFfv34MGDCAl156iXvuuYfmzZvb5OdwFNaUeHPg\neLHlE+QVeyGlVDNgjNZ6sFKqxHNCCHMKCgpi165dzJo1C09PT5vtd+PGjYUlbu3c4aXnCi9vDvBL\nly4V/gUB8MEHH/C3v/0Nd3d32rRpw4oVK266ubLZ2erslEVA8bHyMidqAZgzZ07h47CwMMLCwmwU\nQQhhS1prbty4UVjgR44cKZw7fPbs2Vy4cIEdO3Zw9epVhg4dWvi4f//+hfOEN2rUiMjISEaMGMGl\nS5fQWqOU4sqVK9StW9fqucPLmisc4OjRo/j7+xd+3/Lly3nuueewWCxERUUBeRNw+fj4EBISQmpq\najW/a7YRFRVVmL8y1pT4SaBVseUW+euK6wEsV3l/EzUEhiulsrXWkaV3VrzEhRCO69ixY9xxxx2F\ny4sXL+aNN97g4MGD1KxZk1dffZUlS5aglOKpp54qfDxt2rTCecLr16+Pr68v9evXJzk5GXd3d6ZO\nnVo4vm3tkXhZc4VfunSJiRMnEh0dDeQddb/wwgv84x//QGvN1q1befrpp3nrrbdo1qwZSikmTpxY\nPW+WjZU+wJ07d26521Y6Fa1Syh1IBu4CTgMxwAStdWI5238CbNBarynjOZmK1saMnIr2u+9unjtF\npoV2Xm+//Xbh3OEPPPAAv/zyC127diU3N5f9+/cXPk5MTKRz585YLBYOHjxYOI/4lClT+PTTT+nY\nsSO9evWy6RCNs7utqWi11hal1JPAt+SdzfKR1jpRKfVo3tP6/dLfctuJhRAO54knngAoMZ5doPgR\ne8Hj69evs2jRohLziM+YMcNecV2GVWPiWuvNQFCpde+Vs+10G+QSDiApCRYvLlru1g3y79YlRKVk\nHnH7kMvuRblu3Mj7KpCRYVwWIUTZ5LJ7IYQwMbvfY5M5dns5IYRwDnOQGyU7q+o6O+Xo0W1s2/Yv\nAgJG0KXLZGrW9KVPHxg2rGgbOTtFCPuo6OwUGU4RN9m9+wO++up+goJGc+LEThYvbsuaNZM5f/6U\n0dGEEKXIB5uikMWSzTffzOTw4S1Mn74dP79AevZ8nIyM82zd+k/+/e+JjB//fZnzVgghjCFH4qLQ\njz++RFpaIjNm7MTPL7BwvY+PH0OHvoHWmtdff93AhEKI0qTEBQBZWZfZs+d9wsM/pEaN+jc97+bm\nzrPPfs6CBQvYtWuXAQmFEGWREhcA7N79Hu3bD6N+/dblbtO4cSvefPNNJk6cWDijnBDCWFLigpyc\n6/zyy2L69n2u0m0ffPBB7rzzzsKpPoUQxpISF8TFfUHjxiE0bRpq1fYLFixg6dKlnD9/upqTCSEq\nIyXu4rTOZceO1+jXr/Kj8AINGzZk0qRJrFu3pBqTCSGsISXu4pKTN+DpWYvWrQff0vfNnDmTjRs/\n4Pr1q9WUTAhhDSlxF/fzzwvo2/fZMu9xWJG2bdvSvfvd7NnzQTUlE0JYQ0rchaWlnSM1dR8dO46t\n0vdHRDzLzp0LsViybZxMCGEtKXEXtnXrt7RuPRh3d68qfX9g4B34+QUSH7/cxsmEENaSEndhP/64\nmfbth1W+YQX69n2OHTteRSY2E8IYUuIuKjc3l6iob267xNu1u5fcXAvHj0fbKJkQ4lZIibuoPXv2\n0KCBX4VXaFpDKUVo6BT27v3UNsGEELdEStxFbdq0iSFDhttkX126TCIxcTVZWZk22Z8QwnpS4i5q\n8+bNDB58e0MpBerWbUGzZj3YunW9TfYnhLCelLgLunDhAnFxcfTpY7tb14eGTuHrr2VIRQh7kxJ3\nQd999x0DBgygRo0aNttnhw5jiIvbyenTMp+KEPYkJe6CNm3axPDhthkPL+DlVYvBg8fyxRdf2HS/\nQoiKSYm7GK01mzdvtnmJA4wcOYVPP/1UzhkXwo6kxF1MXFwctWrVol27djbfd7duA7h69Sp79+61\n+b6FEGWTEncxP//8MwMGDKiWfbu5ufGHP/xBhlSEsCMpcRcTExNDr169qm3/ERERrFq1SoZUhLAT\nD6MDCPuKiYnh8ccfr7b9h4SE4O3tTWxs7C3/snjvPcjOnxDRwwMee6waAgrhZORI3IWkp6dz+PBh\nunTpUm2voZQiIiKClStX3vL3nj8PaWlFX0KIykmJu5A9e/YQEhKCl1fVpp61VkGJy5CKENVPStyF\nxMbG0rNnz2p/nS5duuDl5cWuXbuq/bWEcHVS4i6kuj/ULHA7QypCiFsjJe5CYmJi7HIkDjKkIoS9\nSIm7iHPnznHx4kUCAwPt8nqhoaF4eHiwe/duu7yeEK5KStxFxMbG0qNHD9zc7PNPLkMqQtiHlLiL\nsOdQSgG58EeI6icl7iKqcvHN7eratStaa+Li4uz6ukK4EilxF6C1NuRIXCnFfffdx5o1a+z6ukK4\nEqtKXCk1TCmVpJRKUUrNKuP5cKXUPqXU/5RSMUqpfraPKqrq6NGjeHp60rx5c7u/9tixY6XEhahG\nlZa4UsoN+A8wFOgETFBKdSi12Xda61CtdTfgj8CHNk8qqqzg/HCllN1f+8477+TcuXP8+uuvdn9t\nIVyBNUfivYCDWuujWutsYDkwuvgGWuuMYou1gVzbRRS3y4ihlAJubm6MGTOGtWvXGvL6Qjg7a0q8\nOXC82PKJ/HUlKKXGKKUSgQ3AdNvEE7awb98+unXrRnY2JCQUfR0+bNvXOX265P6vXs1bL0MqQlQf\nm01Fq7VeB6xTSvUH5gP3lLXdnDlzCh+HhYURFhZmqwiiHPHx8XTu3JnMTPjqq+p7ndjYvK8CEyZA\nUFDev3NycjKnTp2iWbNm1RdACCcRFRVFVFSUVdtaU+IngVbFllvkryuT1nq7UqqtUspXa32h9PPF\nS1xUv7S0NDIyMmjZsmXhkbG9eXl5MXLkSNatW8cTTzxhTAghTKT0Ae7cuXPL3daa4ZRYoL1Syl8p\n5QWMByKLb6CUalfscXfAq6wCF/ZXcBRuxIeaxcmphkJUj0pLXGttAZ4EvgUOAMu11olKqUeVUo/k\nb3a/UipeKbUHeBN4oNoSi1sSFxdHSEiI0TG49957iY2N5fz580ZHEcKpWDUmrrXeDASVWvdescev\nAq/aNpqwhfj4eIcocR8fH+6++242bNjA1KlTjY4jhNOQKzadnKMciYMMqQhRHaTEnZjWunBM3BGM\nHDmSrVu3ctWoT1iFcEJS4k7s2LFj1K5dGz8/P6OjAFCvXj369evHxo0bjY4ihNOQEndijnQUXkCG\nVISwLSlxJ+ZI4+EFRo8ezTfffENWVpbRUYRwClLiTsxRzkwprlGjRnTr1o1vv/3W6ChCOAUpcScW\nFxfncMMpIEMqQtiSlLiTys7OJiUlheDgYKOj3GTs2LFs2LCB7Oxso6MIYXpS4k7q4MGDtGjRAh8f\nH6Oj3KRFixYEBARYPcGPEKJ8NpvFUDgWRxwPL+6+++5j9erV3HNPmZNd3rKjR2HTpqLloCAYPNgm\nu7a5M2dg3bqi5datYdgww+IIk5MjcSflqOPhBcaNG8eaNWvIycmxyf5u3Mgrx4Kvy5dtsttqkZ1d\nMuvFi0YnEmYmJe6kHP1IvG3btvj7+7N161ajowhhalLiTsrRj8QBHnjgAb6qzrtUCOECpMSdUGZm\nJidPniQgIMDoKBWKiIiw6ZCKEK5IStwJJScn065dOzw8HPtz69atW9O2bVt+/PFHo6MIYVpS4k4o\nKSmJDh06GB3DKhERETKkIsRtkBJ3QomJiXTs2NHoGFaJiIhg7dq1cuGPEFUkJe6EzHQk7u/vT/v2\n7WVIRYgqkhJ3QmY6Egc5S0WI2yEl7mQsFgsHDx4kMDDQ6ChWGzduHGvXriUn54bRUYQwHSlxJ3Ps\n2FEaNWpE7dq1jY5itVatWtGpUydSUjZVvrEQogQpcSeTlGSuoZQCkydPZt++pUbHEMJ0pMSdTEqK\neT7ULC4iIoKDB78lK+uS0VGEMBXHvhpE3LKkpET69OkBgNaQklL0XGamQaGs0KBBA9q3v5uEhNV0\n7/7Hm57PzobDh4uWa9QAf387Bizl8OG8TAUCAsBNDomEAaTEnUxKShJTp04GIDcXli0zONAtCA2d\nzI4dS8os8WvXSv4szZrBI4/YMVwp69eXnCnx73+XEhfGkP/snIjWmuRkc46JAwQFjSA1dT+XLx8z\nOooQpiEl7kQyMtLIzc2lcePGRkepEg8Pbzp2HEdc3JdGRxHCNKTEnUhaWhJBQR1RShkdpcq6dJnM\n/v2fo7U2OooQpiAl7kTS0hIJDDTfmSnFtWrVjxs3rnH69F6jowhhClLiTiQtLYkOHcw5Hl5AKTdC\nQx/if//7xOgoQpiClLgTcYYjcYBu3f7I/v1fkunI50QK4SCkxJ2IMxyJA9Sv70+zZj1YvXq10VGE\ncHhS4k4iOzuD9PQz+Pu3NjqKTfTo8TAffPCB0TGEcHhS4k7i/PkUGjRw/FuyWSsoaBTJyckkJycb\nHUUIhyYl7iTS0pJo2ND84+EFPDy8mDJlCh9++KHRUYRwaFLiTiKvxM0/Hl7cjBkz+Oyzz7hxQ+YZ\nF6I8UuJOwtmOxAECAgIIDg5m/fr1RkcRwmFJiTuJvBIPMjqGzT388MO8++67RscQwmE5x6dgLk7r\nXM6fT8HPr3pLPC4Ojh8vWi4+i19ZvvkGfvqpaDk8HJo0KVpevhyuXi1aLuuG9+PGjeOvf/0rBw7E\nASFVym20K1dgxYqi5evXjcsinI9VJa6UGgYsIu/I/SOt9b9LPT8RmJW/eBV4XGsdZ8ugonyXLx+n\nZk1fvL3rVOvrXLuW92WtCxdKLpcu6dRUuHix4n14eXnxxBNP8O67i2jX7iPrX9yB5OTAyZNGpxDO\nqtLhFKWUG/AfYCjQCZiglCo9+HoYGKi1DgXmA3KCrx0543h4cY899hgbNqzh2rWzRkcRwuFYMybe\nCziotT6qtc4GlgOji2+gtd6ptS7443on0Ny2MUVF0tKSqn0oxUgNGzZkzJgIYmPfMTqKEA7HmhJv\nDhQbCeUEFZf0DEBuW25H588nO/WROMBjjz3Drl3vkJOTZXQUIRyKTT/YVEoNBqYB/cvbZs6cOYWP\nw8LCCAsLs2UEl5Q3Z8pYo2NUqw4dgmnatCtxccvo1m2a0XGEqFZRUVFERUVZta01JX4SaFVsuUX+\nuhKUUl2A94FhWutyP64qXuLCNpx9TLxAnz5/ZsuWZ+nadSpg3htfCFGZ0ge4c+fOLXdba4ZTYoH2\nSil/pZQXMB6ILL6BUqoVsBr4g9b6UBUyiyrKzLzC9etXqFvX+T+GaNfuXkDz668yWidEgUpLXGtt\nAZ4EvgUOAMu11olKqUeVUgX3G/8H4Au8rZT6n1IqptoSixLOnk3Gzy+QvJOInJtSioED/8HWrXPl\n9m1C5LNqTFxrvRkIKrXuvWKPHwYetm00YY3UVNcYSikQHDyOrVv/SXz8JmCE0XGEMJzzH745ubNn\nnfv0wtKUcmPQoNl8/fUcORoXAilx03O1I3GA4OD7uXEjk40bNxodRQjDSYmb3Nmzzn+OeGlKuTFy\n5GzmzJGjcSGkxE0sJyeHtLRD+PkFGB3F7rp1u4+srCy+/vpro6MIYSiZxdDEjhw5Qp06TfD09DE6\nilVOnCg5g19ZsxYW0BoOFTtZNT295PNubm7861//YtasWQwbNgzwrPC1L14sOSGXnx/Ur2999soc\nPgzu7nmPlYK2bW2379LS0krOINm4MdS5hbnPjhwBi6VouW3bvMzCnKTETSwpKYnGjc0zlLJ5s/Xb\nWizw+ecVbzNq1CjefPNN3nnnHYYPf7rCbRMSYMuWouV774W+fa3PU5kvvyx67O0NL7xgu32XtmsX\n7NxZtHz//RByC7P0rlpV8pfiSy9JiZuZlLiJJScn06SJeUrc1pRSLFy4kMGDB9O790SgodGRhLA7\nGRM3saSkJJcucYDOnTszfvx4liyZbXQUIQwhJW5iecMprnOOeHnmzp3LN9+sIjVV7kMiXI+UuElp\nrTlw4ABNmwYbHcVwvr6+PPnkbDZvflpOORQuR0rcpM6ePYtSijp1GhsdxSE8+OAj3Lhxjd2736t8\nYyGciJS4SSUkJBAcHIyS0woA8PDwYOzYz/jxx39w4YJMpClch5S4SRWUuCjSsGEH+vd/gfXrp5Gb\nm2t0HCHsQkrcpKTEy9a7958A+PrrxQYnEcI+pMRNSkq8bG5u7owe/QmrV79MfHy80XGEqHZS4iYl\nJV4+X992TJv2BmPHjuXixXLvFCiEU5ASN6G0tDSysrJo1qyZ0VEcVljYHxg5ciTjx4/HUnyiECGc\njJS4CSUmJsqZKVZ47bXXsFgsvFCdE5kIYTCZO8WEZCjFOh4eHqxYsYKePXtSo0YX3NwmGx1JCJuT\nEjeZ//4XVq5MwM8vmI8/LjklqStJS4OPPy5azswsezs/Pz8iIyMZNOhu7r67Fh07jgUgJgaSkoq2\nGzAAAhx0WvZ9+2D37qJlGeYXxUmJm0xqKvz2WwJNmgzj2DGj0xjnxg2s/vk7d+7Ma69t5Omnh+Pm\n5k5QUDiXLsGlS0XbXLtWPTlt4coV639W4XpkTNyEzp1LoFEjGU65FUFB3Zk48b9s2PAwKSn/NTqO\nEDYjJW4yGRmXyMq6TL16LY2OYjrNmvVgwoQNrF8/jb17/5/RcYSwCRlOMZlTpxJp1KgjSsnv36po\n3rwXU6duZfny0Zw5s5d7730dNzdz/N9Aa016+mmysi6hdS65uRY8PX2wWNoC7kbHcw1Xr0J8PHTu\nfGv3xKtG5vivVxQ6fVqGUm5Xo0YdmTHjF1avnsDSpcMYN245jnhXoN9++40ffviB1au3s2/fAdLS\nkvD0rEnNmn4o5YabmzvXr1/h44/TCA0N4Y477mDy5Mn07t3b6OjO6erVvE/ADxyATp1g2zaHKHJl\nz/mXlVIy2bMQwpT6AFsBLwBPT/jpJ+jTxy6vrZRCa13mhSF2PxKXSftvT0jIcDp3foKgoFFlPj9r\nFtSsmffYYoF58+wYzoF07QpjxhQtR0eXvFFycadO7Wb79kf43e/q8/bbbxMUVPndkhYuLP/0ztI3\nSr5wAZYsKXvb69evkJ6+kdTUNXz77bd07tyZ8PBwfv/73xde0LVtG3z/fflZit8o2WKxEBkZyYIF\nC8jIyGDZsmU3/Tyvv37zjZLdZHSucvlH4tf37cM7ODjvaNwByD+dyZw6JcMpttas2R189tkvjBo1\nir59+zJ+/Hh2Fz8x28auXDnJ7t3v8+WXI3njjRZERy9l6NChpKSksH37dp577jk6depUpSty3d3d\nGTt2LNu2beORRx6hf//+fPLJJ3LwZAt16sC2bQwChxlKARkTN5WrV6+Snn6O+vVbGx3F6Xh4ePDM\nM88wffp0PvzwQ8aOHUu7du2YNGkSI0eOpGnTplXe99WrV/n112iOHNnKoUPfcPnyUdq3H06XLpO5\n//4v6dKlLhMm2PCHIe/P78cee4z+/fszfvx4vv/+ez7++GO8vLxs+0Kupk4dfsn/X0chJW4iCQkJ\nNG3aATc3OROhutStW5eZM2fy1FNPsXbtWtasWcOzzz5LQEAAYWFhdOzYkY4dO3LtWgA5ObVxd/dC\nKUVuroXMzAtkZKSRmXmchQsPEB8fz759+0hKSqJhwx74+w9i2LDFtGx5p93OiOncuTOxsbE8+OCD\nTJo0iWXLliH/t3cu8q9pInv37qVVq25Gx3AJnp6ePPDAAzzwwANkZ2ezbds2du7cyQ8//MBbb71F\nQsIhrl+/Rm5uDh4eNbBYblCjRn18fBpSt+7vqFWrM7169WL69Om0aXMH779fw7CfpWbNmnz11VeM\nHj2aadOmERLyKTKS6jykxE1k7969tGwZanQMl+Pp6cmQIUMYMmRI4bqCDzZzcy3k5GTh4VGj8C+k\nsj7YNFqNGjVYu3YtI0aM4ODBxxg69D2ZBdNJyK9jE8kr8a5GxxDFuLm54+VVyxRDXD4+PmzYsIFT\np/axffsrRscRNiJH4iZhsViIj49n0qRQ0tLK3+7YsbwjQQBXvldwejocOVK0fKtHw6mp5c+MCJCT\nU6VYhqtTpw4PPbSaRYt64O8/kFat+nPkSNEphm5u0KqVoRHFLZISN4lDhw7RqFEjfHzqVbjdsmV2\nCuTgfv0176uqtmy5ve93ZPXrtyA8/CNWr57Io4/+j88+8yt8rmbNvGsNhHnIcIpJ7N27l9BQGQ8X\nthEY+Hs6dXqA9eunyjnkJiclbhJ79+6la1cZDxe2c9ddL3Pt2ll++WWx0VHEbZASNwkpcWFr7u5e\n3Hffl/z003wuXTpidBxRRVLiJrFv3z4pcWFzvr7t6NPnz2za9LTRUUQVWVXiSqlhSqkkpVSKUuqm\njz2UUkFKqR1KqSyl1Ezbx3RtZ8+eJSMjg1Zy2oCoBn37/pXz51NISlpvdBRRBZWWuMq7+8B/gKFA\nJ2CCUqpDqc3OA08Br9k8oWDfvn2EhobKxRmiWnh4ePP737/N5s1Pc/16euXfIByKNUfivYCDWuuj\nWutsYDkwuvgGWus0rfVuwKRnzzo2GQ8X1a1NmyH4+w/k++//aXQUcYusKfHmwPFiyyfy1wk7kfFw\nYQ/33PM6u3d/QnJystFRxC2QDzZNQM4RF/ZQu3YTBgz4Cy+++KLRUcQtsOaKzZNA8U/UWuSvq5I5\nc+YUPg4LCyMsLKyqu3IJmZmZHDp0iOBguRGEqH79+v2Jd98NYOfOnfSx063HxM2ioqKIioqyaltr\nSjwWaK+U8gdOA+OBiqawr/DTt+IlLip34MABAgMD8S6YEEWIauTpWZO5c+cya9YsoqKi5MN0g5Q+\nwJ07d26521Y6nKK1tgBPAt8CB4DlWutEpdSjSqlHAJRSTZRSx4E/Ay8qpY4ppWrf1k8hAPlQU9jf\nlClTOHfuHJs2bTI6irCCVRNgaa03A0Gl1r1X7HEq0NK20QTA7t276dZNbgQh7MfDw4NXXnmF559/\nnqFDh+J5MgnQAAALF0lEQVTu7vjT7LoymcXQwUVHRzN9+nSjYzi96GjYt69o+cwZ47Lcqp9+gj17\nyn9+1Cjw9b21fYaHh/Paa6/xxRdf8NBDDxWuT02FzZuLtvP3B/lYy1hS4g7s8uXLHD58WIZT7ODc\nubwvM6os+40bt75PpRTz589nxowZTJw4EQ+PvKrIyoLffivaroZxd50T+eQUQwe2c+dOevTogaen\np9FRhAsKCwujZcuWfP7550ZHERWQEndg0dHR9O3b1+gYwoXNnTuXefPmkZ2dbXQUUQ4pcQe2Y8cO\n+vXrZ3QM4cIGDhxImzZt+Oyzz4yOIsohJe6gcnJyiImJ4c477zQ6inBxc+fOZf78+XI07qCkxB1U\nXFwcLVq0wPdWTysQwsb69+9P+/bt+fTTT42OIsogJe6goqOjZShFOIyCo/EbVTnVRVQrKXEHtWPH\nDvlQUziMvn37EhAQwJo1MjbuaKTEHZQciQtHM3v2bN5662UsFhkbdyRS4g7oxIkTZGRkEBAQYHQU\nIQr179+fli3bsH+/nDfuSKTEHVDBUIrMICcczZ/+NJtt2/4lR+MORErcAclQinBUvXsPpF69VsTF\nfWF0FJFPStwByYeawpENGjSbn36aT26u3FLXEcgEWAa7fBmuXClaPn/+HCkpB2natCfHj9+8/fXr\n9ssmqiY3lxL/dsX/fa2Rng4XLxYtX7p0e3nOnIHi1+lYLOVva7GUzO7hAb/7XcltWrcOo27dFuzf\nv5ROnaaWeK509tLq14c6dazPLionJW6w3bvzphItsHfvRpo3v4ulS+VOPmaVnQ0ffVT1709Ohg0b\nbJdn3Trrt71xo2T2hg3hySdv3m7w4HmsW/cQERETAa/C9YmJ8N//lr//ESOgVy/r84jKyXCKg0lJ\n2UBgYLjRMYSokL//APz8gvjhh9v4bSVsQkrcgeTkXOfw4S0EBIwwOooQlRoyZD7r1v2LzMxMo6O4\nNClxB3LkSBSNG3emVq1GRkcRolLNmvWgXbuevPPOO0ZHcWlS4g4kOTlShlKEqURE/JNXX32V9PR0\no6O4LClxB6G1JiVlA0FBo4yOIoTVWrUKYciQISxatMjoKC5LStxBpKbuw93di4YNOxodRYhbMm/e\nPBYuXMjp06eNjuKSpMQdRHLyBgIDR8ml9sJ02rVrx4wZM3jxxReNjuKSpMQdREpKJEFBMh4uzOnF\nF19k06ZNJCbuNjqKy5ESdwBXr57iwoVDtGrV3+goQlRJ3bp1mTdvHm+88Qxaa6PjuBQpcQewf/9S\ngoJG4e7uaXQUIaps2rRpZGamk5Cw0ugoLkVK3GAWSw6xsW/Rq9dTRkcR4ra4u7szc+Yitmx5lhs3\n5JRDe5ESN9j27euoW7cFzZr1MDqKELete/dBtGkzhC1bZhkdxWVIiRtszZrF9O79J6NjCGEzQ4cu\nJCUlksOHvzc6ikuQWQyr2dGjsH172c/99tseTpw4wpgxY+0bSjiUEyfgi2L3WLh82bgspV2+XDJb\n6WlSjh0r+fylS1CjRn1GjfqAyMjpPP54HN7ede0T1kUpe36SrJTSrvbJdXw8rFpV9nPr1k2hYcOO\n9O//vH1DCWEHkZEPo5Ri1Kj3C9c5w1S0Sim7n4GT/5plXkQiwykGSU9PJTk5ku7dHzY6ihDVYujQ\nBRw69A3JyZFGR3FqUuIGiY19m+DgCHx8/IyOIkS18Pauy7hxK4iM/COpqfuNjuO0pMQNcP58CrGx\nb9Gv33NGRxGiWrVo0YdhwxazbFk46empRsdxSlLidpaba2HduqkMGvQSvr7tjY4jRLULCZlIaOgU\nVqwYw/XrWUbHcTpS4nb2889v4OHhTa9eZdy4UAgnFRY2m3r1WjF79mSuy92+bUpK3I7OnUtgx45X\nCQ//GKXkrReuQyk3xoz5FK01I0aM4LIjnUdpctIkdnLjxjXWrZvC4MHzadCgjdFxhLA7D48avPzy\nVwQHBzNw4EBOnjxpdCSnICVuB+npqXz66WAaNw7hjjseMTqOEIZxd3dnyZIlTJo0ib59+/Ldd98Z\nHcn0rCpxpdQwpVSSUipFKVXmpAhKqSVKqYNKqb1Kqa62jWlehw8n8tFHfQgMHEl4+Edy0wfh8pRS\nPPfcc7zzzjs8/PDDTJw4kTNnzhgdy7QqLXGVN3j7H2Ao0AmYoJTqUGqb4UA7rXUA8CjwbjVktZmo\nqKhqf42cnBy++OILpk0LIyxsLoMGveSQBX7kSJTRERyGvBdF7PFejBgxgvj4ePz9/QkJCWH+/Pmc\nOnWq2l/3VtmjL26HNUfivYCDWuujWutsYDkwutQ2o4HPALTWvwD1lFJNbJrUhqrzHyUjI4O33nqL\nwMBA3n33XRYuXE1o6EPV9nq3S4qriLwXRez1XtSqVYtXXnmFrVu3cvz4cTp16sSYMWNYv349ly5d\nskuGyjh6iVszAVZz4Hix5RPkFXtF25zMX+fUZ/dfvXqVU6dOcfjwYX7++Wd+/vlnYmJiCAsLY+nS\npfTt25f4eDh82OikQji24OBg3nvvPRYsWMDy5ctZsmQJkydPpm3btvTv35/g4GDatGlDmzZtaN68\nOXXq1HHIv2yNYPdZDEeNGmXvl7xJcnIyu3fn3Quw+EQ2WusSX7m5uVgsFnJycrBYLGRmZpKRkUFG\nRgbnz59Ha02zZs3w9/end+/ePPPMM/Tp0wc/v6JL6WvWhN/9zu4/otVq13bsfPYk70WR6novfHwq\ne93azJgxgxkzZpCdnc2ePXuIjo4mLi6OyMhIfvvtN06dOkVWVhb16tWjfv361KhRA29vb7y9vfH0\n9MTd3R13d3fc3Nxwc3NDKVX4VVzx5Yp+IRTviwKO0GMFKp3FUCnVB5ijtR6Wv/w8oLXW/y62zbvA\nj1rrFfnLScAgrXVqqX251hSGQghhI+XNYmjNkXgs0F4p5Q+cBsYDE0ptEwn8H2BFfulfKl3gFYUQ\nQghRNZWWuNbaopR6EviWvA9CP9JaJyqlHs17Wr+vtd6olBqhlPoVuAZMq97YQgghwM43hRBCCGFb\nLn3FplLqL0qpXKWUr9FZjKKUelUplZh/kdZqpZTL3UvLmovZXIFSqoVS6gel1AGlVJxS6mmjMxlN\nKeWmlNqjlHLYO1u4bIkrpVoA9wBHjc5isG+BTlrrrsBB4AWD89iVNRezuZAcYKbWuhNwJ/B/XPi9\nKPAnIMHoEBVx2RIHFgLPGh3CaFrr77TWufmLO4EWRuYxgDUXs7kErfUZrfXe/MfpQCJ513u4pPwD\nvRHAh0ZnqYhLlrhSKhw4rrWOMzqLg5kObDI6hJ2VdTGbyxZXAaVUa6Ar8IuxSQxVcKDn0B8c2v1i\nH3tRSm0Bil/6r8j7x/g78DfyhlKKP+e0KngvXtRab8jf5kUgW2v9pQERhQNRStUGVgF/yj8idzlK\nqd8DqVrrvUqpMBy4I5y2xLXW95S1XinVGWgN7FN5l2m1AHYrpXpprc/aMaLdlPdeFFBKTSXvz8Yh\ndgnkWE4CrYott8hf55KUUh7kFfjnWuv1RucxUD8gXCk1AqgJ1FFKfaa1driJkFz+FEOl1G9Ad631\nRaOzGEEpNQxYAAzUWp83Oo+9KaXcgWTgLvIuZosBJmitEw0NZhCl1GdAmtZ6ptFZHIVSahDwF611\nuNFZyuKSY+KlaBz4TyU7eBOoDWzJP5XqbaMD2ZPW2gIUXMx2AFjuwgXeD5gEDFFK/S//v4dhRucS\nFXP5I3EhhDAzORIXQggTkxIXQggTkxIXQggTkxIXQggTkxIXQggTkxIXQggTkxIXQggTkxIXQggT\n+/+ldWdIsKG4UwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f309c6c9510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xsample, ysample, accept, c, normalize = sampling(nsamples=1000)\n",
    "plot_samples(xsample, ysample, accept, c, normalize, trace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n"
     ]
    }
   ],
   "source": [
    "# make an animation\n",
    "nsample = 1000\n",
    "xsample, ysample, accept, c, normalize = sampling(nsamples=nsample)\n",
    "divisor = 10\n",
    "\n",
    "for i in range(nsample):\n",
    "    plot_samples(xsample[0:i], \n",
    "                ysample[0:i], \n",
    "                accept[0:i], \n",
    "                c, \n",
    "                normalize,\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": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkz\nODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2Nj\nY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQED\nEQH/xAAbAAEAAgMBAQAAAAAAAAAAAAAAAQUDBAYHAv/EAEwQAAEDAQIGDAkLAwQDAQEAAAABAgME\nBRESFVVylNEGEyExNDVBVHGBscEUIjM2UYKSk7IWMkNhZHORocLh4iNS0iRCRGJTdPBjJf/EABkB\nAQADAQEAAAAAAAAAAAAAAAABAgQDBf/EACoRAQABAwMEAQQCAwEAAAAAAAABAgMRMUFRBBIzcTIU\nITTwE4FCYZEi/9oADAMBAAIRAxEAPwDz8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAHuGIrHyVQ6OzUMRWPkqh0dmosABX4isfJVDo7NQx\nFY+SqHR2aiwAFdiKx8lUOjs1E4isfJVDo7NRruhkq7cq4lq6iKOKKNWtifgpeuFf2IZ8Vfb673wE\n4isfJVDo7NQxFY+SqHR2aiMVfb673wxV9vrvfATiKx8lUOjs1EYisfJVDo7NQxV9vrvfGKpsx7I0\nVlfXX4TU8ryKqXgZsRWPkqh0dmoYisfJVDo7NRGKvt9d74Yq+313vgJxFY+SqHR2ahiKx8lUOjs1\nEYq+313vhir7fXe+AnEVj5KodHZqGIrHyVQ6OzURir7fXe+GKvt9d74CcRWPkqh0dmoYisfJVDo7\nNRGKvt9d74Yq+313vgJxFY+SqHR2ahiKx8lUOjs1EYq+313vhir7fXe+AYisfJVDo7NQxFY+SqHR\n2ajDUWbIx0OBXVyo6REd/V5LlK6niqJLXdTOtCs2tHOTyvoL00TVmY2QuMRWPkqh0dmojEVj5Kod\nHZqIxV9vrvfFZbNPPR7VtNoVqYV998t5FFE11dsJWmIrHyVQ6OzUTiKx8lUOjs1GKCzFfBG5a+tv\nc1FX+sZMVfb673xUTiKx8lUOjs1DEVj5KodHZqIxV9vrvfDFX2+u98BOIrHyVQ6OzUMRWPkqh0dm\nojFX2+u98MVfb673wE4isfJVDo7NRGIrHyVQ6OzUMVfb673xhls2RtTCxtdXYDsLCXbfq3AM+IrH\nyVQ6OzUMRWPkqh0dmojFS8/rvfDFS5QrvfATiKx8lUOjs1DEVj5KodHZqIxUuUK73wxUuUK73wE4\nisfJVDo7NQxFY+SqHR2aiMVLlCu98MVLlCu98BOIrHyVQ6OzUMRWPkqh0dmojFS5QrvfDFS5Qrvf\nATiKx8lUOjs1DEVj5KodHZqIxUuUK73wxUuUK73wDEVj5KodHZqJxFY+SqHR2ajAtmyeGJH4fXYG\n1337by3mbFS5QrvfATiKx8lUOjs1DEVj5KodHZqIxUuUK73wxUuUK73wE4isfJVDo7NQxFY+SqHR\n2aiMVLlCu98MVLlCu98BOIrHyVQ6OzUMRWPkqh0dmojFS5QrvfDFS8/rvfATiKx8lUOjs1DEVj5K\nodHZqIxUvP673wxUvP673wE4isfJVDo7NRGIrHyVQ6OzUMVLlCu98YW2Y9amRi19dgI1qou28u7e\nBnxFY+SqHR2ahiKx8lUOjs1EYqXKFd74YqXKFd74CcRWPkqh0dmoYisfJVDo7NRGKlyhXe+GKlyh\nXe+AnEVj5KodHZqGIrHyVQ6OzUadp0clHQvnir6xXMVt2FLem+hdgSAAAAArKXzhtD7mH9ZZlZS+\ncNofcw/rLMAAABBIAgkAAAAAAAAAAAAIOepPOJ+e46E56k84n57jRZ+NXpEt2stKWntFlO1rVa67\ndXf3TX2SfQdfcYrU48i9XtMuyT6Dr7jrbpimujHCNlxTcGizE7DKYqXg0WYnYZTHOqwACAAAAgkA\nAAAAAAAAAAAAAGOWWOFuFK9rEVbr1W4mOVkqKsbkdcty3cimGtp/CWRsvRESRrnX8qIt935GqrKq\njbM6NrZHzSq/cRd7XciJ1AWQK5KmtdhokVy7l17F5XXdlyny6rq2J47blcqI3xF5XKnZugWZDntY\nnjuRvStxWNkq9sWRYla5ytYt6LciXKt93pvVELB8ayMb4+Cv1NRe0D6bJG9bmva5fqW80orRWaV2\n1Na+NrntW53jeLuKt3SlxtRwqx16yK76sFE7ENKOzXRPjwHNRkU0k7fS5XYW4v1XuX8gIktSRmEx\nIGpLtbZGI59yXOW65dzcN2llfNA179rwl30jdhInWacNJWXvfK6FsjlRyuYqqrrl3t1NxLtzrNmi\npnU+3K9yK6WRZFRN5t/In4AbRBIAAAAAAK63uKJulvxIWBX29xRN0t+JCwAkgFfVV1Q2qlhpYWPW\nCJJJMN11999zU+vcX8gLAGKmqGVNJFUsX+nIxHoq+hUvNOzrSdW1tTFtSMijax8b7917Vv3buTeA\nil84bQ+5h/WWZWUvnDaH3MP6yyAkAAAAAAAAAAAAAAAAAgAc9SecT89x0Jz1J5xPz3GizpV6RJan\nHcXq9pl2SfQdfcYrU47i9XtMuyT6Dr7jvT87fpGy4puDRZidhlMVNwaLMTsMphnVYABAAAACABII\nJAAAAAAAAAAACFVE31uCKipubqGpaNOtVEyNGoqbY1XX8jb938jG6qjo5dojhRGIq3o3oVy3J/8A\nb4G+Q5jXKiuai4K3pem8po4wcmEqxt8ViPW5196ehNzdXcUh1oPa5jdqaquS+5H3rvonf+QFgDVp\nHvqKLDeqK5963ehF3k/AhtK5HIuDFv8A/bWAqa+Ome9qtc7ambZIrf8Aa30/kv4G057Wxq9V8VEv\nv+orK+imfJVrC1HeEwJFff8ANVL91fq8b8jfjw2tcxY/FYiI1b08dLgNfGlPc7CwmqkaSNRUTx2r\nuJd1mxBMsuFfG+NWrdc7l3L9z075UyUU9a2V01MscqqxUVXpvNcioxLl3t/8SwoYZIn1Dnpgskkw\no4/7UuRO1FXrA3AAAAAAAgCvt7iibpb8SFgV9vcUzdLfiQsABWVVNWMrJ5qRsb0qIkYuG7BwHJfc\nu9upu731FoAK9lA5LPbZyrdTtgbGkjXeMqpuLudBjobNmpbUnqH1L5I3xMY1HYPJf6ET0loAKeKB\nk2yGuw1el0MPzXq3+/0KW7URqIicnpKZKd0+yGuunliuhh8mqJfuvNvF7+fVXtJqIyiZnhvA0cXv\n59Ve0moYvfz6q9pNQzPCMzw3gaOL38+qvaTUMXv59Ve0moZngzPDeBo4vfz6q9pNQxe/n1V7Sahm\neDM8N4Gji9/Pqr2k1DF7+fVXtJqGZ4Mzw3gaOL38+qvaTUMXv59Ve0moZngzPDeBo4vfz6q9pNQx\ne/n1V7SahmeDM8N4hVREVV3kNLF7+fVXtJqPmSz37W7/AF1VvL/uTUMzwZnh942of/OnsrqKanqY\nWW0+dz7olc5UdcVG0rf5aT8UMNQ9Ibm7ZK+R3zWNVL3Ho02uyJ+0/f0p3Twvq2eKpteKSF2E29qX\n9Zs7JPoOvuOboZFfWwxSungkc5Lr1RUXoUurepHRbTfVTvvv+cqbn5FczFyiIjROZxouqGtp5mRx\nRyI6RGJely+g3DkrFpnSV2ClRMzxF3Wqmov8Xv59Ve0moy3qeyrGFomeG8DRxe/n1V7Sahi9/Pqr\n2k1HLM8GZ4bwNHF7+fVXtJqGL38+qvaTUMzwZnhvGuytp31CwNkRZEW7BuUw4vfz6q9pNRSU1K5b\ncezwmdFwneMipf2F6Ke6J+2hmeHUg0cXv59Ve0moYvfz6q9pNRTMmZ4bwNHF7+fVXtJqGL38+qva\nTUMzwZnhvA0cXv59Ve0moYvfz6q9pNQzPBmeG8DRxe/n1V7Sahi9/Pqr2k1DM8GZ4bwNHF7+fVXt\nJqGL38+qvaTUMzwZnhvC5DRxe/n1V7Sahi9/Pqr2k1DM8GZ4b1yehD4SJiSrIieMqI3f5P8A5TUx\ne/n1V7Sahi9/Pqr2k1DM8GZ4byXIDRxe/n1V7Sahi9/Pqr2k1DM8GZ4bwNHF7+fVXtJqGL38+qva\nTUMzwZnhvA0cXv59Ve0moYvfz6q9pNQzPBmeG8DRxe/n1V7Sahi9/Pqr2k1DM8GZ4bwNHF7+fVXt\nJqGL38+qvaTUMzwZnhvA0cXv59Ve0moYvfz6q9pNQzPBmeGvb1OzF00l8mFe1fKOu+cnJfcWxSWz\nROjs2R61lQ9EVviuVLl8ZPqLsJhIAJSAACspfOG0PuYf1lkVtL5w2h9zD+sswIBIAgEgCASAIBIA\ngEgCASQAPmTyb+hT6PmXyT+hRGo4lVuvVd41KJNtV9U5N2VfE+pib2s+rQVdoSJN+Z6R/jv/AJXm\nw1Ea1GtS5E3EQ9rWfTmbR4Q+NjXYD8NFY9N9rr9xSxtGqdV0tOsrcGaNzo5W+hyXX9XL1mpScKhz\n07TZ2RxJS2pDKi3MqkwXJyYaby9abnUhwuYi9TKY0ZLA4x9RTpjmbA4x9RTpjL1XkWp0ASDMlAJA\nEHPUnnC/OcdCc9SecL85xos6VekS6EAkzpQCQBAJAEAkAQCQBAJAEAkAQCQBAJAEAkAQCQBAJAEA\nkAV1vcUTdLfiQsCvt7iibpb8SFgBIAAAACspfOG0PuYf1lkVtL5w2h9zD+sswAAAAAAAAAAAAAAQ\nSAIPmTybuhT6PmTybuhRA8/k/qWkxn+2JivXpXcT8rzaNWjXDlqZfTKrU6G7nbebR7dGmXOWWl4V\nDnp2lhsxjWSmiwPntvczpS5UK+k4XDnp2lvsk+g6+4z3YzephMaNLY1Ik1XHI3efEq9h1Jx+xT+n\na09Ov+y9zc12723lpUVM7bdbE2VyR4TfFv3N4z3aZuVz/qExovQAZFgAAQc9SecL85x0Jz1J5wvz\nnGizpV6RLoSSCTOkAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABXW/xRN0t+JCwK+3uKJulvxIWAEgg\nq6qWqmrqiGnn2ltPE16eKi4blv3Fv5NxPxAtQalHULXWZDUNXa3TxNelyX4N6XmpZ1VI+0Zqfb5J\no2MRyrMzAcjr+TcS9APul84bQ+5h/WWRV072M2QV+E5G3ww7655aIqKl6ASAAAAAAAAAAAAAAACD\nHUvSOmle5bkaxVX8DKVeyJypY08bb8Ka6JLv+y3d4HJ0LcGjjvS5XJhL0ru95nCJgpcnJuA9ymMR\nhzZaThUOenaW+yT6Dr7iopOFw56dpb7JPoOvuM9zzUpjRUWZI2mt2mmXelRYF691PzS7rLKq84m5\n7ewrqanWp25jE/qJGr485qoqfmhtMqG1drU9QzelwHfiiFKoxdq9JjR1IAPOWAABBz1J5wvznHQn\nPUnnC/OcaLOlXpEuhJIJM6QAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAIAr7f4om6W/EhYFbb0jMVTN\nwm4V7dy//shZADQq7PklqJJoKhYHSxpHJ4t96Jfcqehd1SwIA1o6V0ESRU8u1xsiSONuDfg3cv18\nhjgopUq/CqmdJZEYrGo1mCiIq3ryr6EN0AVcEbJNkNfhsa66GG69L/7y0RERLk3itpfOG0PuYf1l\nmAAAAAAAAAAAAAAAAAKq2L5KqggTeWV0juhrV71QtSqlVJbdcnN6VV63r/FRGo5td8DlB7rky0nC\n4c9O0t9kn0HX3FRScLhz07S32SfQdfcZrnnpW2atgcY+oprU0fg2yKSl5I50czNdu9t6GzYHGPqK\na1ry+CbKoqm69q3Mcn13Xp3nO5n+accJjR2AKT5Qt5uvtFy1b2ovpMVduqj5Qtl9AHy5bmqvoQoJ\nOepPOF+e4srOtJK572pGrMFL9+8raTzhfnuNNqmae+J4RLoSSCTMkAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAgCSCQBV29FHiuZ+A3Cvbu3bvzkLMr7f4om6W/EhYASAAAAArKXzhr/uYf1lmVlL5w1/3MP6\nyzAAAAAAAAAAAAAAABAElVRptldak/8A3SJF+prdaqWhV2Je6yXzLvzySSdSuW78rhA5td8Bd8Hu\nw5MtLwuHPTtLfZJ9B19xUUvC4c9O0t9kn0HX3Ga556VtmrYHGPqKamyiN0lTPgfPbgvb0oiKbdgc\nY+opjtrjOXq7CJjN+Y/0bK6GRssTJG7rXIiod1H5NvQee0P9NZab/wAT9zNXdTvTqPQo/Jt6EOPV\nTmKUwkiTybug+j5k8m7oMSyg2N+WmzUIpPOJ+c4nY35efNQik84X5zjfX5LnpV0JJBJgWAAAAAAA\nAAAAAAAAAAAAAAAAAAABAJAAAAV1v8UTdLfiQsCvt/iibpb8SFgBIAAAACspfOGv+5h/WWZWUvnD\naH3MP6yzAAAAAAAAAAAAAABAPl8jI0ve5Gp6VW4DDaEyU9BUTKtyMjc6/qMdBD4PZEEP9kKJ+Rr2\n+5JLHfG1UclQ9kW5yo5yIv5Xm9NLFGxzXSNauDvKpMajjQF3we45MtJwuHPTtLfZJ9B19xUUnC4c\n9O0t9kn0HX3Ga556VtmrYHGPqKY7a4zl6E7DJYHGPqKY7a4zl6uwR559Gynm/o10MqfNk/pv7U70\n6z0GPybehDgquJZqZ7G7jt9q+hU3UO2s6pbWUEFQzcR7EW70LyoZ+rjEwmlsnzJ5N3QfR8yeTd0G\nNZQbG/LTZqEUnnC/OcTsb8tNmofNO5rNkEjnORqYTt1VN9fkr9K8OiJPhsjHMwmuarfSi7hDJY5L\n8B7XXb9y3mBZkBBIAEHwyaOR2CyRrlTkRQMgIAEg+HyMjS972tRfStxLXI5EVFvReVAPoEEgAAAA\nAAAAAAAAAAAAAABXW/xRN0t+JCwK+3+KJulvxIWAEkAp6qPwu0qyOWR7WwQNWNGuVuCq3+Nucu4n\n4AXAvNKge6usWmfPejpoGufgqrVvVEv3U3jXsxrfDp30mH4GjUberlVHPRd1Uv8AwA+6XzhtD7mH\n9ZZlPFOyHZBXYaPW+GH5rFd/f6ELdq3oipy+kCQAAAAAAAAAAAAEFTsi4HHn9xbFTsi4HHn9x1se\nSETowIquhsaBP98quXoa1y9txh2Q8NZmJ2mazk22voPRBSvXrc5E/SYdkPDmZnedrHn/AOk6KoAH\npubLScLhz07S32SfQdfcVFJwqHPTtLfZJ9B19xlueelbZq2Bxj6imO2uM5ersMlgcY+opjtrjOXq\n7BH5E+jZonR2b/pK+oorrmSJ4RF1/OT8d3rOc5TpLVa6Kngr4/n0io913LHd4yfhu9Ry6zZNK0Ik\n8m7oDHI9jXtW9rkvRU5UEnk3dBgWUOxvy02ahXWhw+fPUsNjflps1CvtDh8+ep6lHnqU2XNlcSS+\nv2GHY1/yPV7zNZXEkvQ/sMWxr/ker3mer41+0rwkgkxrIXeOesHjGXNXtOhXeOesHjGbNXtNFrx1\nol0IAM6VPsj4PDn9xYWfwCDMQr9kfB4c7uLCz+AQZiHerw0/2jdsEkEnBIAAAAAAAAAAAAAAAAAQ\nBX2/xRN0t+JCwKm3qhmLpo7n4V7U8m675yct1xbADUqrOp6uTbJMNrlbgOVj1bhN9C3b6G4ANaSi\nika9q4bWvj2pUa9URG/V6Ok+aOgiokRsL5cFEwUa6RXIifUim2AKyl84bQ+5h/WWRW0vnDaH3MP6\nyzAAAAAAAAAAAAAAIKnZFwOPP7i2KnZFwOPP7jrY8kInRg2O/wBWoqpeSNkcKdSK5fiQw7IeHMzO\n82dirP8A+W+Xllme78FuTsNbZDw5mYnadrHn/wConRVAA9NRlpOFRZ6dpb7JPoOvuKik4XDnp2lv\nsk+g6+4y3PPSts1bA4y9RTHbXGcvV2GSwOMfUUx21xnL1dgj8j+jZonatajoUa5L0VtyocVynbR+\nTb0Ices2TSrrIVaZ81myLu063xfXEvzfw3U6iyf5N3QV1rRviWO0IG3y01+E1P8AfGvzk7+o3mys\nmpkljcjmPbhNVOVDCso9jflps1CvtDh8+epYbG/LTZqFfaHD589T1KPPV6UnRc2VxJL0P7DFsa/5\nHq95lsriSX1+wxbGv+R6veZ6vjX7SvQQDGsLvHPWDxjLmr2nQrvHPWDxjLmr2mi1460Tq6EAGdKn\n2R8Hizu4sLP4BBmIV+yTg8Wf3FhZ/AIMxDRX4af7Ru2CSAZ0pAAAAAAAAAAAAAAAAIJAFdb3FE3S\n34kLAr7e4om6W/EhYASAAAAArKXzhtD7mH9ZZlZS+cNofcw/rLMAAAAAAAAAAAAAAgp9kzsCga9d\n5rr/AMlLgoNmbrrGVE33OwU69w6WpxXCJfdlvkpdiNNIzckSFHdaru9pU1NVLVyI+ZUVyJduJcXN\n2DsXRqckaJ+ZQG7paY+875VqAAbFUscrHte3fat6GeqrZqvB25UXB3rkuNcFZpiZyllpqiSll2yJ\nUR1126l5FRO+omWWRUVy79xjA7Yz3bhynbR+Tb0IcTynbR+Tb0IYet2WpVPh06214NhJtWFddd9R\nMK4tqnUTuDT3up1/tXfVnenWaqecvr9xcWhSxVlJJDMi4KpeiotytVN5U+s4XqYjtxwmFRsb8tNm\noV9ocPnz1M+xOoclVNT1Cok6MS70SJ6UMFocPnz1Nluc36vSs6Jhr6iCnWBjkRi33pd6Sz2N/wDI\n9XvKMvNjf/I9XvJ6imItzgjVekAHlrqm2q6eklibCqIjkVVvS809jy310irvqxe0+9knl4M1e0x7\nHeGvzO83xTEdPM/uqu7owAYFlPsj4PDn9xWxWtVxRtjY9qNalyeKWWyPg8Od3HPnp9PRTVajuhSd\nXUUFVLPZj55FRXphXLd6DHYtbPWbbtzkXBuuuS4+LJ4kl9fsMWxv6fqM1VFOK/tpKV6ADKsAAAAA\nAAAAAAAAAAArre4om6W/EhYFfb3FE3S34kLACQAAAAFZS+cNofcw/rLIraXzhtD7mH9ZZACSCQAA\nAAAAAAAAAg5/Zfc6kpI/7qhq/gir3HQHM7J3q+0aOFPmxsdI7pXcTvOlqM1xBLfd5s+onac8dC7z\nZ9RO0549Dpv8valQADUqAAAAAB2zPJt6EOJO2Z5NvQh5/W7L0qBvnKuf3F/J5N3QUCecvr9xfv8A\nJu6Djf8A8fSYcxZNFHWula9XMexEdHI1bnMX0oVlXUSU9dLHW3YWGqJMiXNd0+hS82N+XmzUK60m\no6uqGuRFRXreipvmqnP89WFdmuioqXoXuxv/AJHq95zPgSMW+nlkh/6t3W/gvcXFgNtRu37TLSOT\ncv2xjk9PoUt1Ez/HOYI1dUCtwLafuLPQxp/1ic5fzcgdZc06/wCrtGpe3lZEqRN/Ld/M8tdWbKam\nCGeDbJWNXBXcv3fwNSwq2Z1Y/wAFo5Zb2fOeuA38938jYtqgpaSaLaIGMVWre669V6V3zLsd4a/M\n7zfET9P+8q7rLCtl6XpHRR/Ur3O7kIWptOn3aiijmYm+tM+93srd2lkDAs5+2K2CtpI3QPvVr7nN\nVLnNW7eVN9CnLbZVRteyKeG6OpRbkkTlT0L6UKOmqNuaqObgSMW57F5F1Hp9LP8A4iFKnT2TxLL6\n/YYtjf0/UZbJ4kl9fsMWxv6fqOFXxue08L0AGNYAAAAAAAAAAAAAAABXW9xRN0t+JCwK+3uKJulv\nxIWAAxzVMECtSaVkauW5qOciXmQpanaMZ2j4Xg4PgjbsPewPGvu6+4C6VyIiqq3InKfEU0UzcKKR\nsjb7r2rehqWWqusmljqVRZfB2LI12/vbt5isVY2rWsjwURKl1yN6gJpfOG0PuYf1lkc/PjH5Q1mL\n1iu2mLC2xt/993KhutitlURVrKRF9Hg6/wCQFoCs2m2ee0mju/yG02zz2k0d3+QFmCs2m2ee0mju\n/wAhtNs89pNHd/kBZgrNptnntJo7v8htNs89pNHd/kBZgrNptnntJo7v8htNs89pNHd/kBZgrNpt\nnntJo7v8htNs89pNHd/kBZHL2rdNW1k3/jeyH8Gq5fiQtXQ2zgr/AK2k3ubu/wAjlbLW0JLMqZ66\nVsjZatytVOVUvRV6Nw62PJCJ0Xa2jBibwXxtswbt7c3ypAPVooijON1MgAOiAAAAAAOkZblGjERd\nsvRP7Tmwcrlmm58kxOFpTzNqLebKy/Bc69L+g6KTybuhTlbI4zh6e46qTybuhTB1UdtURHC0KLY3\n5abNQrrQ4fPnqWOxvy02ahXWhw+fPU1UeepE6Ncs7HroaLbduwvHuuuS8rAd66IrjtlETh0uPaP/\nAPT2T7itmlmlZGzDwnLcl7TlzYs/jCnz0M1XS0RTMpysdknl4c1e0x7HeGvzO8ybJPLw5q9pj2O8\nMfmd5WPxv3lO7owAeesp9kfB4c7uOVqo3RuSqiS97E8dqf72+g6rZHweHO7jnz0+njNqFJ1X9jPb\nJYL3sW9rkcqL1Hxsa+n6jTsOR1JTzUb0uhma+SBfr/3N7zc2N/T9RnnPbczylegAyLAAAAAAAAAA\nAAAAAQBX29xRN0t+JCwKq3pJMXTN2lcG9vj4Sf3IWoEmGamgnVqzQskVi3tVzb7jMAMaxMVyuVjV\nVzcFVu309B8Q0lPTuV0MEcbl3L2tRDOAKyl84bQ+5h/WWZWUvnDaH3MP6yzAAAAAAAAAAAAAANS0\n51prNqZkS9zI3K1PSt24n4nM0sC0uxuKBVvWOpkaqry3OUvrbVZG0lKm/PUMRc1vjL+TbusqZ+K5\nP/em+Nx1seSETorwAey5gAAAAAAAAAA3LI4yh6e46qTybug5WyOMoenuOqk8m7oU8zrPnC9Ki2N+\nWmzUK60OHz56lhsb8vNmoV9ocPnz1NNHnq9I2a4ANSobFn8Pgz0Nc2LP4wp89Clz4SmFjsk8vDmr\n2mPY9w1+Z3mTZJ5eHNXtMex7hr8zvMkfjfvK27owAeesp9knB4c/uOfOg2ScHhz+4589XpfHClWq\n2jpFqtjr8BcGaJzpInehyJ373WfWxOVJ4ZJE3MJEvT0LyobNk8Sy+v2FZsZkSktKanVLo6lMNi8i\nPTfTr3/xMtecV+0urBBJkWAAAAAAAAAABAAAkAAV1vcUTdLfiQsCvt7iibpb8SFgBIAAAACspfOG\n0PuYf1lmVlL5w2h9zD+sswAAAAAAAABBJCgL0F5x9RNKlRKiSPuwl/3L6TorGcrrMiVy3ru7q9Kn\ne5Ym3TFWUROWN67fshjYnzaaBXLnOW5PyRSpn4sl/wDem+NxbWPfLHU1jt+olcrcxPFb2X9ZydK5\nypOiuVU8Il3L/wDspPTU91cTwTozAA9ZzAAAAAAAAAABuWTxlD09x1T/ACbug4pFVq3tVUVOVDob\nDc59nyq5yr4y76/UYOro/wA16Zauxvy02ahX2hw+fPUwMe5nzXK3oUhVVVVVW9VNNNvFya+VdgAH\nZAbFn8YU+ehrmxZ/GFPnoUufCUwsdknl4c1e0x7HuGvzO8ybJPLw5q9pj2PcNfmd5kj8b95W3dHe\nLzlbVlkbaMyJI5Ev3kX6jU26X/yP9pTnT0k1RE5O5ebI+Dw5/cUB9Oe9/wA5znXelT5Ntqj+OntV\nmcuhsriWX1uwrrPpXVVHUJEt08atkiX0OS+78d7rNFJHtbgo9yIvIil1sb+n6jNdt9lFc8piVrQ1\nTaykjnaiphJutXfavKi9CmwVUa4vtd0S7kFauGz0NkRN1OtN3qUtTz1wAAAAAAAAAAAAAAAFdb3F\nE3S34kLAr7e4om6W/EhYASAAAAArKXzhtD7mH9ZZlZS+cNofcw/rLMAAAAAAAAAQpJCgcXU8Jlz1\n7S4ZO6HY61sXlplWKLOcqp+W/wBRT1PCZc9e0s7DRa2aJ9y7TRo5Gqu86RV3buhNzrU9DqfHCtOq\n7poG01JHAz5sbEanUhwtHvT/APsSfEp36/NU4Cj3p/v5PiU5dJ8ypsAA9NQAAAAAAAAAAA6GweLp\nc5ew546GweLpc5ewy9X41qdXPAA1KgAAGxZ/GFPnoa5sWfxhT56FLnwlMLHZJ5eHNXtMex7hr8zv\nMmyTy8OavaY9j3DX5neZI/G/eU7ta1uM5+nuNM3LW4zn6e40zVb+EIkAB0QF5sb+n6ijLzY39P1G\nfqfFK0arO0aNK2kfDhKx+45j032OTdRfxPmzaxaymvkbgTsXAlj/ALXJv6zcKuujfR1OMadiuTBw\naiNu+5qbzk+tOw8ldZkmOGVk8TZYnI9j0va5N5UMgAAAAAAAAAAAAABXW9xRN0t+JCwK+3uKJulv\nxIWAA1Ku0aekk2uTDVyNw3YDFdgt9K3byG2U9VKlHaVY+WN7mzwNSNWtV2EqYXi7nLup+IFq6WNs\nKzK9NrRuFhX7l3pMFJaENW9WR7Y12CjkSRitVW+lLzQWCT5O4tuctSyjaipctyrdddfvch9U0iVl\nrU88LXpHDTuY/CarfGcrbk3fRgqBlpfOG0PuYf1lkV0tBVYwlqqWqZFtrGtc18WF82/609JPg9q8\n/g0f+QFgCv8AB7V5/Bo/8h4PavP4NH/kBYHzI9sbcJy3Jeifitxo+D2rz+DR/wCR8S0dpyswXV8N\n16LuU/oW/wDuAswV/g9q8/g0f+Q8HtXn8Gj/AMgLALvFa6ltZyblowt6Kf8Ac0Z7DtSoeizW5I5n\nLGkSNavTcqKBQ2hUOkq5oKVUWTCXCfyR7vb9R1mx6JkFjQRxpc1t/XuqVjdi8zEujrYmN9DYLk7T\nfpqC0qaFsUdfDgt3r4P5Gi7dprpiN0RC1XeU4Ck+n/8AYk+JTr/B7VX/AJ8Gj/yKqPYtUR4eDaDf\nGer1vh5VW9eUdPci3VmSYyrQWnyaqcoM9z+4+TVTlBnuf3Nn1dtXtVTnI26/lW5CSzdsXqHK2+0G\n+Kt6f0f3J+TVTlBnuf3I+rtnaqwWnyaqcoM9z+4+TVTlBnuf3J+rtnaqwWnyaqcoM9z+4+TVTlBn\nuf3H1ds7VWC0+TVTlBnuf3HyaqcoM9z+4+rtnaqzobB4ulzl7DS+TVVlBnuf3Nulsu0KSF0UdoRY\nKret8H7nC/foroxCYjCgIVyI5Gqu6u8WvyaqcoM9z+58rsXqFe1y2g29t939H9zv9XbR2q0Fp8mq\nnKDPc/uPk1U5QZ7n9x9XbO1VmxZ/GFPnobnyaqcoM9z+59w7HquGVsjbQZhNW9L4P3K1dVRNMwdr\n62SeXhzV7THse4a/M7zYrLHrq1zXS2hFe1Lkug/c+aSxa6jkWSK0I71S7dg/c4Rep/h7N04++Wha\n3Gc/T3GmXFRYFZUTOlfaDMJ2/dB+5j+TVTlBnuf3O1HVURTEI7VWC0+TVTlBnuf3HyaqcoM9z+5f\n6u2dsqrCTDwb9268vdjf0/Ua3yXqcPDxg2+67yP7m3R2TX0WFtVoReNv3wX95xvdRRXRNMJiF2QV\n/g9q8/g0f+Q8HtXn8Gj/AMjCsxPp57NldNQsWWmet8lMi7rV9LNRuUdfT1jb4ZEVyfOYu45v1Km+\nhg8HtXn8Gj/yNSqsWqq3I+aqp9sTekbT4L06FR14F2Chgsa2IPm29I5v9r4Ud+a7puNpbWTftKFe\nmn/cCyBX+D2rz+DR/wCQ8HtXn8Gj/wAgLA+Ue1ZHMRfGaiKqdJo+D2rz+DR/5HwlHaaSukSvhwnI\niL/p/Rf/ANvrAswV/g9q8/g0f+Q8HtXn8Gj/AMgLAFf4PavP4NH/AJDwe1efwaP/ACAW9xRN0t+J\nCwKqps+0aqFYZq+JY3Kl6NguXcW/0lqABIAgEgCASAIBIA+JJGRRukkXBa1L1X0GrHalFJI1jJ0V\nzluRLl3TcFwFHPaEktrUixVCMptudErUVP6i4K3qv1IqXfiXhozWRRzVEMywxtdE/DS5ieMtypu/\nibwAEgCASAIBIAgxzzx08aySuwWJvqZSANSK0qSokSKGdHPXeREU1WbZDa0MEVTNP4quqEeqKjU5\nF+pb+TpLVUvT0GhQ2a+ikVyVksjXOVz2va3xlX0rdeBYAACASAIBIAgEgCCkrXTSvtKVtRJGtHHf\nEjVubhI3CvVOUuzRqrLiqpZHrJKxJWoyVjV3JETkX9gNqmkWamikVLlexHKnShkDURqIibiJuISB\nAJAEAkAQCQBAUkgCpqttgrqdsFVNJNLLe6JyorUjv8bcu3ETk+u4tivisx8NbJUtrZb5X4Tmq1q3\np/bfdfcWIEAkAQCQBAJAEAkAQq3JepSQV8lRbVO5s6JSysejIkVPGuuucvTu3fUXapely7xpYpo0\nrIalkMbHxIqIjWIm/du/kBugkAQCQBAJAFfj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj\n5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM\n1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8\nq0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7\nhj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZr\nPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5V\nodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1j\nHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0\nOkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj\n2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPD\nwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5VodIZrPDwB7hj2x8q0OkM1jHtj5Vod\nIZrPDwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAB//9k=\n",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/hObkFth7kMs\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7f13723c5e10>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "YouTubeVideo(\"hObkFth7kMs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\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
}