LuxXx/MetropolisHastings.ipynb

## MetropolisHastings.ipynb
{
 "cells": [
  {
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# may be unnormalized\n",
    "def target(x):\n",
    "    if x < 0:\n",
    "        return 0\n",
    "    return np.exp(-x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# may be unnormalized\n",
    "def target(x):\n",
    "    return 1/(np.sqrt(2*np.pi)) * np.exp(-x**2 / 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# may be unnormalized\n",
    "def target(x):\n",
    "    return np.exp(-x**2 / 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# may be unnormalized\n",
    "def target(x):\n",
    "    return np.sin(x)**2 * np.sin(2*x)**2 * np.exp(-x**2 / 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.zeros(10000)\n",
    "x[0] = 3\n",
    "for i in range(1, len(x)):\n",
    "    curr_x = x[i-1]\n",
    "    prop_x = curr_x + np.random.normal(0, 1)\n",
    "    a = target(prop_x) / target(curr_x)\n",
    "    if np.random.uniform() < a:\n",
    "        x[i] = prop_x\n",
    "    else:\n",
    "        x[i] = curr_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5840424474631826"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = quad(target, -np.inf, np.inf)[0]\n",
    "N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x194f8310>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYXHWV8PHvqeruJN3pTtJLNhLoAAEJiDBGFlEGFccwanAfcJlxRmV8X2HcxpkgyoyoMwz4uiC4oDI6iDKIiIFEA8ouAZNAgISAJDEknXSS3qv37qo67x/3VnWl00t19626S5/P8/CkqvrWvccyOXX6/JYrqooxxphoifkdgDHGGO9ZcjfGmAiy5G6MMRFkyd0YYyLIkrsxxkSQJXdjjIkgS+7GGBNBltyNMSaCLLkbY0wElfh14draWq2vr/fr8sYYE0pbtmxpVtW68Y7LK7mLyCrgW0Ac+KGqXjvs598A3uA+LQfmq+rcsc5ZX1/P5s2b87m8McYYl4i8nM9x4yZ3EYkDNwFvBhqATSKyVlWfzxyjqp/OOf4K4MwJR2yMMcYz+fTczwJ2qupuVR0AbgcuHuP4S4GfexGcMcaYycknuR8D7Mt53uC+dhQROQ5YBjww9dCMMcZMVj7JXUZ4bbR9gi8B7lTV1IgnErlMRDaLyOampqZ8YzTGGDNB+ST3BmBpzvMlwIFRjr2EMVoyqnqzqq5U1ZV1deMO9hpjjJmkfJL7JmC5iCwTkTKcBL52+EEicjIwD9jobYjGGGMmatzkrqpJ4HJgA7ADuENVt4vINSKyOufQS4Hb1W7tZIwxvstrnruqrgfWD3vt6mHP/927sIwxxkyFbT9gjDER5Nv2A8YUUv2addnHe659q4+RGOMPq9yNMSaCLLkbY0wEWXKPioYt8MBXoC2vPYWire1lPlNyB2fKS35HYoxvLLlHQfNL8JO3wyPXw0/fBYN9fkfkn4EeuPWd/FPJ3fy87CsskcN+R2SMLyy5h1D9mnXZ/wB49OuAwsXfgZadsPWnvsbnq623Qesu/nnwH0kT41Mld/kdkTG+sOQecqes+SV9W3/Brb2vhTM/AAteCU9P4+T+1P/AoldxZ+ovuTv1Wi6KPQkD3X5HZUzRWXIPufNjzzJTBlmfPtt54ZXvgQNPQ2K07X8irPMQHHwWTn0XAGvT51Eh/bDz9z4HZkzxWXIPufNjz5HQWfwx/QrnhRPf5Py560H/gvLLnx92/jz+AgC2pE+iW2fAnx/xLSRj/GLJPeT+IvYSW9MnkiLuvLDgNCivhZf/4G9gftj9MMyaBwtPB2CQEjalXwF7HvU5MGOKz5J7iFXQy0myj6d1OeAOtF65ngc6lzitmelm70Y47jyIDf21fiq9HJpehP5OHwMzpvgsuYfYq2K7iIs6CSzHc3o8NL0wvQYS+zuhdRcsOuOIl7dpPaBwcJsvYRnjF0vuIXa67Abg6fQJR7z+bHoZaBoan/UjLH9kkvei0494+bn0MudB4zNFDsgYf1lyD7HlsQYatZoEs494/dlMsp9OrZmDzzl/LnzlES83MQ9mL7DkbqYdS+4hdpI08FL66HuVNzEXZlU7rZnp4uAzzkBy5aKjfvRgx0Kef/oPR+wUaUzUWXIPKSHNiXKAP+mSkQ+oO9nZlmC6OPwCLFgBcvT93HfqMSyTRoS0D4EZ4w9L7iG1VJqYJQOjJ/fa5dD8YnGD8osqtLwENctH/PFuXcQsGWAxLUUOzBj/WHIPqeXSAMDOEdoyANSeDD0t0D0NElpPC/R1OF9oI9idXgzA8bHGYkZljK8suYdUvRwCYJcuHvmAupOdP5v/VKSIfNSy0/mz5sQRf7xbnT788WLJ3UwfltxDaqkcJqHldFAx8gGZKnY6JPfM2MIoyb2JOSS0nONlGu63Y6YtS+4hdawcZq/OB44eQASgaglIHNr3FjUuX7TshFgpzD12lAOE3bqQZXKwqGEZ46e8kruIrBKRF0Vkp4isGeWY94nI8yKyXUR+5m2YZrih5D6KeAnMOQbap8GdmVp3w7x6iMVHPaRB57NEmooXkzE+Gze5i0gcuAm4CFgBXCoiK4Ydsxy4EjhPVU8FPlWAWI1LSLNUmtindaMeU79mHY+3zGbLM1uLGJlPOvaNUbU7GrSOY6QZ0jYd0kwP+VTuZwE7VXW3qg4AtwMXDzvmY8BNqtoGoKp2b7MCmk87M2SQfWNV7sA+nc/S6VCttu/NI7nXMkOS0HWoSEEZ4698kvsxwL6c5w3ua7lOAk4SkT+IyBMismqkE4nIZSKyWUQ2NzVNg6RTIMe69wUdsy0D7NM65ks7DPYWIyx/DHQ7UyHnLh3zsAatdR507BvzOGOiIp/kPtKInQ57XgIsBy4ALgV+KCJzj3qT6s2qulJVV9bVjd5SMGPLN7k3ZNo27RFOaJn/bXOPG/Owoc9iGgwwG0N+yb0ByC2LlgDD55Q1AL9W1UFV/TPwIk6yNwVwbOwwaRX2j9FzB4Z68lEeVM0k6zljV+77M5V7lD8LY3Lkk9w3ActFZJmIlAGXAGuHHXM38AYAEanFadPs9jJQM+QYaeYQ8xikZMzjsj35tj2FD8ovHW5yH6fn3stMWrTSKnczbYyb3FU1CVwObAB2AHeo6nYRuUZEVruHbQBaROR54EHgc6o6Dda9+2MhLTRq9bjHNTGHfi2NdrXavhfiZc62vuNo0DpL7mbaGLv0c6nqemD9sNeuznmswGfc/0yBLZJWXtCx2xAASowDWs2yjv1FiMon7ftgzpIjbq03mgat5VVRHn8wJoetUA0bVRZKKwe1Jq/DD2oNJCK87L5977j99owDWuvMltHh8wGMiR5L7mHT10GF9OfVlgE4yDzojHByT+x3Kvc8HNR5kOyD3rYCB2WM/yy5h41bhR/MN7lrNSQao7kyM51yFiWNcPelkRzKfGadtseMiT5L7mHjJve8K3ethvSgs9AnarqbnBuBV+WX3A/qPOdBlH+TMcZlyT1sEs7g6IQq95z3RUqnuz97npX7QTKfhe3rbqIvr9kyJkASB0ircIh5eR0+VK02AmcULq4iytzo+sLYFn5YBlQuzOt9h4/4LIyJNqvcwyaxnybmkMzze7kxM6smgpX7AnEHRitHuRvVMAOUwqxqS+5mWrDkHjaJA3n32wGamePctCOCrYj50gYSg4oJ7FNUtTiSn4Uxw1lyD5vOg0PthTykiTltiwjOdV9IG1TMd25Mkq/KRTagaqYFS+5h03WIJj1qw82xRTShLZC2vGfKZFUtsqmQZlqw5B4mKWdKYxNzJva+qkWRbEUskLa8Z8pkVS6CrsPOZ2lMhNlsmTDpbgZ04pX77IXw50cLEpKfFkjrETNlMrNoxlS5CFBn8VOeK1uNCSOr3MPEvUVck06wcp+9APraIdlfgKD8UcYg1dKV90yZrMyXgd1uz0ScJfcw6XLuwDTxyn3+Ee+PgvnS7jzIc457VgQ/C2NGYsk9TLrd5M5Ek7u713mEEtoCWp0HE+y5v/bb2wD41//5vdchGRMoltzDZNJtmUy1Gp1WxNACpolV7i1UAVBHu9chGRMoltzDpOswzJhDP2UTe1+2co9Ocq+VDufBBJN7P2V0aPnQ+42JKEvuYdJ1CGZPYDVmRmYFZ4TaMjWSIK0Cs/Jf0JXRrHMsuZvIs6mQYdJ1OK97hR6lpMzZUyVClXsNCVqpZOXnfzvh9zYxlzpL7ibirHIPk65DQ/3ziZq9IFrJXRK0auWk3tusc6jFkruJNkvuYTLZyh2cL4UItWWqpZOWiQ4su5p0DnViA6om2vJK7iKySkReFJGdIrJmhJ9/WESaRGSr+99HvQ91mhvshf6EVe6uWjqyM18mqknnUCW9MNjncVTGBMe4PXcRiQM3AW8GGoBNIrJWVZ8fduj/qurlBYjRwFDVPdnKvXKBcw5VEPEuLp/USIKW9CTbMpm9eboPw9xjPYzKmODIp3I/C9ipqrtVdQC4Hbi4sGGZo0w1uc9eAMle6O/0LiaflJBkrnRPui3TnHlfhNpUxgyXT3I/BtiX87zBfW24d4vIsyJyp4gs9SQ6MyTTUplKWwYikdDm4XxBTb4t467wjcBnYcxo8knuI/0Or8Oe3wPUq+rpwO+An4x4IpHLRGSziGxuamqaWKTTXSa5V0w2uUdnlWqtJABo0ckl96HKPfyfhTGjySe5NwC5lfgS4Ig7P6hqi6pmthz8AfDqkU6kqjer6kpVXVlXN4nFONNZ12FAoKJ2cu+viE5yr55ics9W/N1WYJjoyie5bwKWi8gyESkDLgHW5h4gIrm7N60GdngXogG47febaNHZ1F913+ROkFml2t3sXVA+qcFN7pNsywxQSrtWWFvGRNq4s2VUNSkilwMbgDhwi6puF5FrgM2quhb4JxFZDSSBVuDDBYx5WqqWxFA7YTLKqwGBnvAn96m2ZcBpzcyNwG8xxowmr+0HVHU9sH7Ya1fnPL4SuNLb0EyuaumkdQrJjFgcymsi0YqolgSDGidB+aTP0UIVJ/a0eBiVMcFiK1RDooYELUxuXndWRW0kkntmXxmdwl/fVq2MRIvKmNFYcg+JaklMrXIHp+/eHf5qtdaDz6JVqyLRojJmNJbcwyCdYi7dtE61co9QW6Z5ism9hUroaYV0yqOojAkWS+5h0NNKTHRKA4iAU7lHoFp1WlQeVO4o9LZ5E5QxAWPJPQzchDz1tkytk8xSgx4E5Z8ar9oyYH13E1mW3MPATUBTbstkFkCFeZbIYB+V0jvltkxzpvKPwG8yxozEknsYuAnIk7YMhLtazfwW40lbhnB/FsaMwZJ7GHR71JYpdyv3MA+qdnvzRdeSuYuTVe4moiy5h4HbRmlj9tTOk6ncw9yW8Si5t2VaXBGYGmrMSCy5h0F3M206mxTxqZ2nIgqVuxP7VGfLJCmBmXOscjeRZck9DHqaJ30z6CPMnAsSD3dy92r8AZw2lfXcTURZcg+D7uYpV6oAxGLuQqYQJ7TuJvq1lC5mTf1cFbVWuZvIsuQeBj0ttHlRuYO7BUGIE1p3s7vHjgf3gS2vtZ67iay8doU0PutupkVPnfTb69esyz7ec0rIq9XuZm9aMgAVNbB/szfnMiZgrHIPunQaelqmPK87K+w7Q3Y3TfrG2Ecpr3VmDunwu0YaE35WuQdY/Zp1zKGLZ2ampj7HPSMSbZljvTlXRS2kk9DXDrPmeXNOYwLCKveAq8nedcijnnt5LfQnINk//rFB1NPsbeUO1nc3kWTJPeCq3fuFetqWgXBW7wPdMNjjbc8dwj0GYcwoLLkHXI10Ah5sPZCR3TwshAnNHSuY8gZqGeUh/qIzZhyW3AOu2uu2THbzsBAOqrrtkyndKDxXmL/ojBmHJfeAq8at3D1ry2SSewj7zJmtB7z6LcYqdxNheSV3EVklIi+KyE4RWTPGce8RERWRld6FOL3VSIJOncUApd6csNztM4excvdou9+s0plQNjvcG6kZM4pxk7uIxIGbgIuAFcClIrJihOMqgX8CnvQ6yOnMuTG2Ry0ZcDbLipWGM7lnK3cPP4+wb8dgzCjyqdzPAnaq6m5VHQBuBy4e4bgvA9cBfR7GN+1V0+ldpQogEt49VbqbobSCXmZ6d86wfhbGjCOf5H4MsC/neYP7WpaInAksVdV7PYzN4LRlPK1UwV2lGsKE1t08NH3RK7YzpImofJL7SDs0Zddri0gM+Abw2XFPJHKZiGwWkc1NTSFsC/igWjq9mwaZUR7SLQi6m4YGhL1SUWs9dxNJ+ST3BmBpzvMlwIGc55XAacBDIrIHOAdYO9KgqqrerKorVXVlXZ3H/0gjSakm4W1bBsK7BUF309AMF69keu62v4yJmHyS+yZguYgsE5Ey4BJgbeaHqtqhqrWqWq+q9cATwGpVte32pmg2vcyQZGHaMmGsVntaClO5p/phoMvb8xrjs3GTu6omgcuBDcAO4A5V3S4i14jI6kIHOJ1Ve706NaO8xklmg73enreQVN22jNeVu811N9GU166QqroeWD/statHOfaCqYdlAGrcfWVavFpun5FdyNQMc5eOfWxQ9CcgNeB9cs+uUm2B6mXentsYH9kK1QDLbD3geeUexhtlZyprr9syVrmbiLLkHmDZtkwhBlQhXH33bHL3ui1T7fxpc91NxNjNOgIss6+M5wOq2S0IQpTQMr9llNcC+z05Zf2adVTQy/aZhOuzMCYPVrkH2DzppE9L6WWGtycO486QPYVpy3Qzk34tDddvMcbkwZJ7gNWQoIUqRl5HNgUzKiFeFq5WROaLyOu2DOIMWFtyNxFjbZkAmyedtHnckqlfsw6AjTNmsyhMrYjuFphRBSUe/xaDM2C9OEyfhTF5sMo9wGqk09sdIXO0aFW4+syFmOPuatXKcP0WY0weLLkHWHW2LeO90CW0Quwr42qlMlxfdMbkwZJ7gBWiLZPRzJyQDai2eL+vjKtVq6CntSDnNsYvltyDKtlPlfR6d0u5YVq1Mly32itgW6ZFq2CgE5L9BTm/MX6w5B5UbiXZ5vXWA65WrYLBbhjoKcj5PZVOu3u5F7AtA9aaMZFiyT2o3H645wuYXM2ZXn4Y+u597aCpwg6oQjg+C2PyZMk9qNx514XquWf3qwlD371Q+8q4hj4LS+4mOiy5B5WbaAo1Wybbyw9D3z279YDHt9hzZdsyNqhqIsSSe1Bleu6FmuceprZMgbYeyMh+0YXhszAmT5bcg6qnmbQK7cwuyOmzfeZQtGUyWw8UJrl3UAESs7aMiRRL7kHV00I7FaQL9H9RF7MgPiMcCS0TY2Z7Xo8pMZhVbZW7iRTbWyaoupu9v0nHEcSZfRKS5N6msznzqvsKd42wfBbG5Mkq96DqaRka6CuUitpwVKvdTQVbzJVVXmsDqiZSLLkHVU9LgSt3nIQWhmq1u7lgs4ayKmrC8UVnTJ4suQdVTwutWpjB1KywtCJ6motTuYfhszAmT5bcg0jVbcsUulqtC0e1WpS2TA30tkE6VdjrGFMkeSV3EVklIi+KyE4RWTPCzz8uIs+JyFYReUxEVngf6jTS1wHpZMHmuGeV18BgDwx0F/Y6U5FOQU9rEb7oagG1vruJjHGTu4jEgZuAi4AVwKUjJO+fqeorVfUM4Drg655HOp24Ww8UvFrN3ks1wNV7TyugNBejcge73Z6JjHwq97OAnaq6W1UHgNuBi3MPUNVEztMKQL0LcRpyk22hdoTMymzEFeTWjLuAqeCDy2H4LIyZgHzmuR8D7Mt53gCcPfwgEfkE8BmgDHijJ9FNV9nKvdBtGTehBblyd5N7s84p7HXC8FkYMwH5VO4ywmtHVeaqepOqngD8K/CFEU8kcpmIbBaRzU1NIVj27he3eixatRrkhJZJ7oXuuWfbMgH+LIyZgHySewOwNOf5EuDAGMffDrxjpB+o6s2qulJVV9bVFWafkEhwK/eiLGKCYO8vU7TK3U3uYdgl05g85JPcNwHLRWSZiJQBlwBrcw8QkeU5T98KvORdiNNQdzOUltPHjMJep2w2lMwMdrXa3QSxEhKUF/Y6JWUwY44NqJrIGLfnrqpJEbkc2ADEgVtUdbuIXANsVtW1wOUiciEwCLQBf1fIoKOsfs06vlb6LOfECpzMAETcxTsBTmhdh6G8Fu0pwpIMW6VqIiSvjcNUdT2wfthrV+c8/qTHcU1r1SSGtuQttIqagLdlmmF2HRQj59oqVRMhtitkAFVLZ+EXMGUEeJVq/Zp1/KrsJTp1VnEuWF4DHfvGP86YELDtBwKomkThN8rKCHi1WksHzRR4MDWjoibQn4UxE2HJPYDmSVcRK3c3uWsw153VSKLwM2Vwfkv47qYOBjqbAvtZGDMRltwDpoxBKqW38AuYMipqIdkbyP1lZtFHufQXfhsGV4tWUSYp6E+Mf7AxAWfJPWDm0QkUYeuBjMz+MgHsu9dKB0DRWlRD95UN3mdhzERZcg+YGnGqxoKvTs3ILrsP3nTIWpzPohhtGchZNGZz3U0EWHIPmHniVO7FmwoZ3FWqmS+6gu8I6cp+oVrlbiLAknvA1LhtmYJvPZAR4N0Qs22ZolXubnIP4GdhzERZcg+Y6mxbpvDJvX7NOlb81xYArr3z0YJfb6Jq3LZMsXru2UFsa8uYCLDkHjA1kiCpMdop8P1TXT3MoFfLqHbbQUFSKx0ktJwBSotyvV5m0qtl1pYxkWDJPWBq6aCFKrRo/9cILVRl+9tBUisdReu3Z7RSaZW7iQRL7gFTK4mi9ZgzWrUy2wIJkhoSxVud6mrVSqvcTSRYcg8YP6rVJp2bHbwMkhpJFG0BU0arVtmAqokES+4BU1PMvVRcTTqHOmkv6jXzUSsdRU/uLVRZW8ZEgiX3gPGjWm1iLrV0QDpV1OuOKZVkHl0+tWUsuZvws+QeJAPdVBRxL5WMwzqXuGiwKtaeFmKiRVudmtGqlTDYDYO9Rb2uMV6z5B4k2ZtBF7stM9d50HWoqNcdk/tZFL3njq1SNdFgyT1I3IRS/AFV98uk05J7dvGYDaqakLPkHiRdh4HiLbfPOEwAK3c3lmxsRZL9MglSi8qYSbDkHiSZtkyRk3v2ekFK7p0HATis84p6WWvLmKiw5B4k3U7lXrRNw1y9zHTuUxqk5N51iC6dSQ8zi3rZ7BddAHfJNGYi8kruIrJKRF4UkZ0ismaEn39GRJ4XkWdF5Pcicpz3oU4D3c0ktJx+yop+6cM6N1jJvfPg0FhAMS/LLCiZGazPwphJGDe5i0gcuAm4CFgBXCoiK4Yd9jSwUlVPB+4ErvM60Gmhu6nog6kZzczJ9vwDoesQhyluS8YhMHt+sD4LYyYhn8r9LGCnqu5W1QHgduDi3ANU9UFV7XGfPgEs8TbMaaLrcNG2tx2uKZCVe3EHU7NmLwjWZ2HMJOST3I8B9uU8b3BfG81HgN9MJahpq7u56IOpGYd1brCmQnYdcmLyw+wFVrmb0MsnucsIr+mIB4p8EFgJXD/Kzy8Tkc0isrmpyQasjtLdVPR53RlNOhcGOmGg25frH6G/Cwa6/EvuFXVWuZvQyye5NwBLc54vAQ4MP0hELgSuAlarav9IJ1LVm1V1paqurKurm0y80ZVKQk8LLUVenZrRlLluECrWzBx3Pyv3nlZIDfpzfWM8kE9y3wQsF5FlIlIGXAKszT1ARM4Evo+T2AOQHUKotxVQX2aIQE4iDULFmpnj7suAKs6AKmpz3U2ojZvcVTUJXA5sAHYAd6jqdhG5RkRWu4ddD8wGfiEiW0Vk7SinM6PJrk71sS0DwUjuXZkFTD5W7hCMz8KYSSrJ5yBVXQ+sH/ba1TmPL/Q4runHXcDk14DqUHIPwC9enQFoy0AwPgtjJslWqAZFthXhT0JrpRIkno3DV10HIV5WtJuEH2X2fDcOq9xNeFlyD4pOf1sRaWJOUgtCcu885FbPI03UKgJL7iYCLLkHRdchmFFFb5H3UjlC1WLoPGoiVPF1HRxqjfihdBbMCNiKXWMmyJJ7UHT6nNDASe6JACT3zkNQudDfGGbPt8rdhJol96DoPOh/Qqs6JhjJPbHficVPsxfYzpAm1Cy5B4XfrQhwKvf+BPQl/IuhL+HEULXYtxDq16zj3t1Jdv15N/Vr1vkWhzFTYck9CFSD0YrIVMudjf7FkPnNYY6/e8816RzqpN3XGIyZCkvuAfDKK++EZC9feaTN30Ay1XJiv38xJBrcWPxtyzTpXKqkl5mMuJOGMYFnyT0A5ouT1H1btON63fdeBOBzt/i4qWemcvexLQNDe+3USoevcRgzWZbcA2C+++t/k08LmDIy9ytdSKt/QXTsBwQqF/kXA0MrdudjrRkTTpbcA2A+wajcByilSatYJP4l99t//wSHdQ71X7jftxgAGrUagIU+fhbGTIUl9wDIVO6H1KddEHMc0mpfE9piaeGAm1j9lEnufn7RGTMVltwDYL6006Mz6GKW36HQqNUsFP8GdhdJK41a49v1MxJU0KtlVrmb0LLkHgALpM1tyfi0l0qOg1rNQmnx7foLpZWDAajcQWjUaqvcTWhZcg+A+dLu226QwzVqNdXSBYO9xb94XweV0suBAFTu4HzRLbDkbkIqr/3co2y0FYh7rn1r0WJYRAvP6AlFu95YDmit86CjAWqXF/fi7jTIYFTu0Eg1Z8sLfodhzKRY5e63dJqF0hqYanWfuve2bdtT/It3OIungvJZHNJqFtAG6bTfoRgzYZbc/dbTzAxJBmIQEWCfunuZ+5Lc9wI5vz34rFGrKZUU9Ni9VE34WHL3W4ez3L4xIK2IJubQr6XQ/nLxL962h34t4aBfN8YeJtse8nM7BmMmadr33H2XCFYrQomxT+s4sc2P5P4yDVqHBqTmyH7hJg7A4jOLcs3cMaBijvuY6AnGv6LpzB1EDEpbBty+ux9tmbY9Q22hADiUm9yNCRmr3Ecwnzb4w7egfR/MfwW88n0ws6owF+tooF9LaaFA55+EfTof2v9Y/Au37WGfriz+dUfRTBWDGqe0kMm9vwue+Tk0vQBzj6OWWprdTcuMmYq8kruIrAK+BcSBH6rqtcN+fj7wTeB04BJVvdPrQIvl4thjfLX0Fri/z7mPZn8HPPoNuOQ2WHyG9xdM7Hd//fd/AVPGPq2Dvg7obYdZRZp/39sOfe3sDVDlrsQ4xDyWFCq5H9gKP7/E2T/f/bv28IwZXDn4UdamzyvMNc20MW5bRkTiwE3ARcAK4FIRWTHssL3Ah4GfeR1gMV0Sf4BvlX2H7VoPl2+BNS/DP2wAEbj1HdC80/uLduwPVEsGGEqwxRxUda8VpLYMuH13d9DbSxdc+UMS37+IhsQg/MN9cOVeuHwL23QZN5TdxHvjD3l+TTO95FO5nwXsVNXdACJyO3Ax8HzmAFXd4/4stBOCz41t56slP+LB1Kv4+OCn6f/ai4Czv/lS+Qy/LvsiB254N+8cuIZB92PzZMArsZ8DHDf183ioIXc65KJXFeeibcFM7g1ax2va93p70tQg3y39FkliXDLwRR479mzn9doT+eDA5/lR6fVcW/ID2P1WOP4Cb69tpo18BlSPAfblPG9wX5swEblMRDaLyOampuDcfLiGDr5d+m1262I+MfhagDdEAAASkklEQVRJ+ik74uf7dAFXDn6M02J7+Nv4fd5dOJWExIHAVe5DC5mKWLm7A7hBasuA+1kkGiA16N1Jn/w+p8T2cuXgx2jIfNauQUr4+OCn2a2L4a7LoNvm2JvJySe5j9QM1slcTFVvVtWVqrqyrq5u/DcUyRdLb6WSHv7v4CfpYeaIx2xIv4ZHUq/kipJfMYcuby6caABNBS6hJaiAmXOKNmOmfs06bvvNQ7TpbDopL8o187VP54Omef0XbvXmZtndLfDQf/JA6gw2pEcePO5hJlcMXuGMQ6z/3NSvaaalfJJ7A7A05/kSIDJzw86PPcM74o/z3dRqXtKxb8r81eQHqKKHj5V48I8csskzaK0IAGpOhJYCjDGM4nhpZJf6e2u9kWQq6yXiUQX95HdhoIv/SL6fsQbRX9Bj4fWfhe13wa4HvLm2mVbySe6bgOUiskxEyoBLgLWFDatIkv18ueS/2ZVexHeTq8c9/EU9lvvSK/lg/HfMom/q18+0ItIBTO61J0HzS0W73AmxA+xKBy+5Z1pUS+XwlM912ppfkHj4Jn6Teg07xykkAHjdp6D6BFj3WUgOTPn6ZnoZN7mrahK4HNgA7ADuUNXtInKNiKwGEJHXiEgD8F7g+yKyvZBBe2bzLRwXO8yXkn97VJ99NDcn38pc6ea98Yenfv22PRAroZFg9dwBZ0fIzgPQ31nwS1XRzXxpZ5f6e9/UkTRqDUmNsVSmPkb07vijVElPXoUEACUz4KLroHU3bL5lytc300te89xVdT2wfthrV+c83oTTrgmPvg54+DoeS53KI+nT837bU3oSW9PH88H470DVmSY5WW17YO6xpHsCuFC49iTnz+Y/wTGvLuiljpdGAGcQMWBSxGnUmqlX7qpcGn+AZ9LH82ye2zs7PX7lp6Wn8rqH/wvOuNQZCzEmDwHMKoVXv2YdN3zln6C3lWuTlzLRBUT/m3oDJ8X2c/Hnb6B+zbrJD7S17YF59ZN7b6HVnuz8WYTWzAniDOEEsecOTt99ypV7wyZeEdvHz1JvmuAbxfk72tvqrJo2Jk/TMrnX0cZH479hbepctunxE37/Palz6dUy/maqC02CnNyrl0GsxKncC+z42AEGND40BTNg9nmR3Lf8mC6dyT2pcyf81m16PJz2Htj4HUg0Ti0OM21Mi71lhlfWXym5i1KSfC35vkmdr4ty1qfP5u3xjVyT/BB9zJj4SXrbobctsMm9/qr7+H1ZHS899AirJlpsTtAJ0sjLupBkQP867tGFzJeHqWCStx4c7IXnf829qXNGnGqb129+b/wCPP9rePhaeLtV8GZ8065yP04O8jfxh/h56o3s1QWTPs+dqfOplF7eFHt6cifIzCEPaHIHp02yXLxfej/cybKXl3RS6+KKIjPQmxkbmLCX7oOBLtamXzv5IKqXwcq/h6duLcw2GCZypl1y/3TJnSSJ8+3kO6Z0nifTp3BY5/K2+MbJnSDT7sj0tgPoeT2OZXIQBroLd5G+BMtih9ieri/cNaYoMxaQGRuYsG13QUUdT6ZPmVog538OSmbCA1+e2nnMtDCtkvsrZC+rYxv579QqmqZ4t580MdalzuYNsa2T+nX9xjvuJakxln/9xSnFUUjb0suIicLBbYW7yCHn3M9rsPbXybVXF5DUGCfEJpHc+7vgTxtgxcWkiE86hvo166j/yia+1fsWeP5u2P/UpM9lpofIJvfMLJbcfuZnS35BF7P4XvJtnlzj3tQ5zJRBLoxtmfB7T5QD7NGF2U3IgihbTR98tnAXaXz2yGsF0ACl7NX5HD+Zyv1Pv4VkL5z6Lk9i+UHyr6G8Bn7/JU/OZ6Irssl9uL+QP/Hm+Ba+l3wbCWZ7cs6ndDkHtJq3xZ+Y8HuXSwM7A9xnBmikmhathMathbvIwWdp0ioOU6R94ydply7mhMn03LfdBZWL4NiJz5IZSRfl8Pp/ht0Pwa4HPTmniaZpktyVfy65gyat4sepVR6eNcb61Nn8ZewZZ/ZLvpIDHCeHAj2I6BCnom4sbOX+fLqeIN2sZCS7dDHLpBHSqfzf1NcBO++HU98JMQ//qb3mIzDnWPjdv0M6tLtsmwKbFsn9jbGneW38eW5MvnPUXR8n697UuZRJCl6YwEKm1l2USJqdAdxLZbjtWg+Hd0Cy3/uTD/ZC0w7nGgG3WxcxQ5IT2ynzhfWQGvCsJZNVMgPe8Hlo3MonvvBvU1tIZyIrUsl9pD57GYN8seRWdqYXc9uEVweOb6uewL50HWz7Zf5vanIGUfPaPMpnT6WXQ3oQ9k98XGFc+7dAOsnm9Enen9tjL6SPdR4cmsDg8va7nAp7SQHuC3v6+2D+Cj5bcgclJL0/vwm9SCX3kXw4/luWxQ5xTfJDBVokI9yTPtfpgXa35PeWg8+R1Bg7A7rcPtcf068ABPb8wfuT73WmkW4JQXJ/UZeS1Fj+LaqeVmer3lPfMbX9h0YTi8Obrub42EE+FL8fGLm4MdNXpJP7fNq4ouRufpc6k0fShbtd3D2pc0FTsOPX+b2hcSsv6ZK8d6L0UwezYcFp8PJj3p/85Y1Qd4pzjYDrp8wZAM935tCOtZBOwmnvLlxQJ63iwdSr+GzJL1hEnoWFmTaCOw9vypSvlt5CKUm+kvxgQa+0Q49lZ3oxTb/+AZfeuWDse6uqwoGneS59WkFj8lT9ebDlJ86e4iUefSElB2Dfk3D63xx5E8cA267H8Yp8K/fn7oSaE6m/oQHY72kcuZX5Evl77i/7F/699Cf84+BnPL2OCbfIVu6rY4/z5vgWvpZ8H3sKvk+4cE/qXM6O7WA+bWMf2tEAPS08O4kNy/zyj4/OgmQv7/nit7076d6NMNAFJ17o3TkLbHt6GXQdHHfzrrPW/JT0nx/jm4dOp9CzgBp0Pt9Mvpu3xDezKvbHgl7LhEskk/siWvhS6U94On0it6QuKso170mfS0yUt443533/ZsBZ/RkWj6dPZUDjXBj3cFXkS/dBvAyWne/dOQtsS3q582Dv2FtOvC3+BDFR1qamsJfMBPwodRHPppdxbekPWIzdUNs4IpfcS0lyY9kNlJLkM4P/h3SR/ifu1sVsS9ezery9ZvY8BqUVbAvB9L+MTsp5PH0aq2KbnLbSVKnCi+vhuPNgRvD77RnbtR5Ky8dN7qvjj7MtXV+0m48kKeGKwSuIk+aGshtt9owBIpfclS+V/JhXx17iXwcv489Fvm3b2tS5nBnbOeoNLurXrOOlJ3/DQ30nBnZ729FsSK+kPnYIGp+Z+skaNju3jivkYGMBJCnh0b7j2fHEb0efkXJ4B2fEdnF36ryixvayLuSqwY+wMvYnrin5sTdfwibUIpXcPxm/i/eXPMB3kqtZlz6n6Ne/K3U+gxqHLT8e8ed1tLM8tp/H0yuKG5gH1qXOpk9L4amfTP1kz/yMPi3ltDtmhm7a3pPpUzhZ9lHHKCuSt/yEAY1zV+r1xQ0MWJt+Ld9Jrub9JQ/AI9cX/fomWEKf3J15vffy7S/8HZ8u/SW/SJ7Pdcm/8SWWZuawIb0Stt7mrL4c5k1uz/qx9CuLHdqUJZjNPalz6dr0M05fc8fkk3J3MzxzO/emz3X2SQmZ36VfTUyUC+MjLOoa7IVnfs6G9Gtopar4wQHXJ9/HL1Ovgwe/Cg981Sr4aSz0yb2CXr5ZehNXlNzNz5NvYE3yY/i5T8ltqQudOyw9c/tRP3t7bCO70osCvb3tWG5JXUQ5/Xy85J7Jn2TjjTDYy3eTb/cusCJ6QZfycno+b4ltPvqHW38Gfe38NPnm4gfmUmJ8bvDj/G/yAnjkOu6++q+dbYfNtJNXcheRVSLyoojsFJE1I/x8hoj8r/vzJ0Wk3utAR/TS/dxb9nneHtvI1wbfy5XJj05pz2wvbEyvgGNeDY993ZnLndGxn3Niz3Nv+hyCvknWaHbocdydPo+PxH/DSTLxyel/deX3GHj0Bu5KnceuwG+aNhphXfocXhd7Dtr3Dr2cHIDHvglLzuJJfYV/4eHca2BN8qN8bfC9vD22Eb7/etj5O19jMsU3bnIXkThwE3ARsAK4VESGN40/ArSp6onAN4D/8jrQrOSAs0DkR2+B296DIrx/4AvcmHonwUiaAhdc6fzDf/yGoZc33ogi/CL1l/6F5oH/HHw/Ccr5Xuk3IDGB/c0Tjfyg9P+RoIIvDxZ2UVmh3Zp8M4rAxpuGXvzDt6BjL1ywhiD8PVRi3Jh6Jx8YvAokBj99t/NvZvuvRmwZmujJZ8rGWcBOVd0NICK3AxcDz+ccczHw7+7jO4EbRURUC9Dwe+Q6Z7BoXj2s+i/ecveiwN3wov5HA9xUehYX/v4/mbH4TGc+9x9v5pep82nQ+X6HNyVNzOXjA5/ix2XXwffPh/P/BZa/GeYe6+x3kpFKQqofOg86VePD11ErnXxg4CrafOpHe6WRGn6Veh3v2/RDOGU1DPbS/8C13J8+m8t/2Od3eEd4Ir0C/s/jsPm/4YnvwC8+TK+W8WT6FC644M2wYAVUHcNffnc7bTqbAUp54T8uPvL/y6DJTSs5j5d93hkHUoQ9//nWwuzpMwW541RjrmL3iIyXf0XkPcAqVf2o+/xDwNmqennOMdvcYxrc57vcY0ZdUbFy5UrdvHmEvuV42vfy4ev+h4fTr0IDPGQwl07uKLuGk2LO0vNd6UW8a+BLodhHJR8ny17+o/RHvDqWM+0zVgqxEiep65H7jG9On8SawY+GYifMfMyhi1+VORt3AbyYXsL7Bq4O9P+/MdK8NradN8We4rWx7Zwcb3T2RBpBUmMkiTOzdKKFk9I36JxTcl4DmFESyybjwdTQ34/S2NBxqbQe8d6YeFMfplVQIB6TobPLWI/dKEZ8PDld/UPrD2ZffD38xd9O6jwiskVVx91qNJ/k/l7gLcOS+1mqekXOMdvdY3KT+1mq2jLsXJcBl7lPTwaG30C0FgK5xM7imhiLa2IsrokLamzFiOs4Va0b76B8vpYbgKU5z5cAw5utmWMaRKQEmAO0Dj+Rqt4M3DzahURkcz7fSMVmcU2MxTUxFtfEBTW2IMWVT19jE7BcRJaJSBlwCbB22DFrgb9zH78HeKAg/XZjjDF5GbdyV9WkiFwObADiwC2qul1ErgE2q+pa4EfArSKyE6div6SQQRtjjBlbXqMlqroeWD/statzHvcB7/UgnlFbNj6zuCbG4poYi2vighpbYOIad0DVGGNM+AR3LqExxphJC1RyF5Evi8izIrJVRO4TkUDcQVpErheRF9zYfiUic/2OKUNE3isi20UkLSK+jtKPt02FX0TkFhE57K7HCAwRWSoiD4rIDvf/w0/6HROAiMwUkT+KyDNuXF/yO6ZcIhIXkadF5F6/Y8kQkT0i8pybuyaxgMd7gUruwPWqerqqngHcC1w93huK5H7gNFU9HfgTcKXP8eTaBrwLeMTPIPLcpsIvPwZW+R3ECJLAZ1X1FOAc4BMB+cz6gTeq6quAM4BVIlL8PbRH90lgh99BjOANqnpGmKZCFo2qJnKeVpBZ2uYzVb1PVTPLy57AmesfCKq6Q1WHLwbzQ3abClUdADLbVPhOVR9hhHUXflPVRlV9yn3ciZOwfN9RTR2ZrSRL3f8C8W9RRJYAbwV+6HcsQReo5A4gIl8VkX3ABwhO5Z7rH4Df+B1EAB0D5G4V2UAAElVYuDupngk86W8kDrf1sRU4DNyvqoGIC/gm8C9AerwDi0yB+0Rki7sS33dFT+4i8jsR2TbCfxcDqOpVqroUuA24fOyzFS8u95ircH6Vvq1YceUbWwCMtPFGIKq9oBOR2cAvgU8N++3VN6qactujS4CzROQ0v2MSkbcBh1V1hDul+O48Vf0LnLbkJ0TE9zu/F307RVW9MM9DfwasA/6tgOFkjReXiPwd8DbgTcVefTuBz8xP+WxTYYYRkVKcxH6bqt7ldzzDqWq7iDyEM2bh94D0ecBqEflrYCZQJSI/VVXf95BW1QPun4dF5Fc4bUpfx8EC1ZYRkeU5T1cDL/gVSy4RWQX8K7BaVXv8jieg8tmmwuQQEcFZ3b1DVb/udzwZIlKXmREmIrOACwnAv0VVvVJVl6hqPc7frweCkNhFpEJEKjOPgb/C/y/CYCV34Fq33fAszgcUiKlhwI1AJXC/O9Xpe34HlCEi7xSRBuBcYJ2IbPAjDnfAObNNxQ7gDlXd7kcsw4nIz4GNwMki0iAiH/E7Jtd5wIeAN7p/r7a6VanfFgEPuv8ON+H03AMz7TCAFgCPicgzwB+Bdar6W59jshWqxhgTRUGr3I0xxnjAkrsxxkSQJXdjjIkgS+7GGBNBltyNMSaCLLkbY0wEWXI3xpgIsuRujDER9P8BtSd1JqHhvXsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(x, bins=100, density=True)\n",
    "f = np.vectorize(target)(np.linspace(np.min(x), np.max(x), 1000)) / N\n",
    "plt.plot(np.linspace(np.min(x), np.max(x), 1000), f)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
	{
	"cells": [
	{
	"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"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 2,
	"metadata": {},
	"outputs": [],
	"source": [
	"# may be unnormalized\n",
	"def target(x):\n",
	" if x < 0:\n",
	" return 0\n",
	" return np.exp(-x)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 3,
	"metadata": {},
	"outputs": [],
	"source": [
	"# may be unnormalized\n",
	"def target(x):\n",
	" return 1/(np.sqrt(2np.pi)) np.exp(-x**2 / 2)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 4,
	"metadata": {},
	"outputs": [],
	"source": [
	"# may be unnormalized\n",
	"def target(x):\n",
	" return np.exp(-x**2 / 2)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 5,
	"metadata": {},
	"outputs": [],
	"source": [
	"# may be unnormalized\n",
	"def target(x):\n",
	" return np.sin(x)*2 np.sin(2x)2 np.exp(-x**2 / 2)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 6,
	"metadata": {},
	"outputs": [],
	"source": [
	"x = np.zeros(10000)\n",
	"x[0] = 3\n",
	"for i in range(1, len(x)):\n",
	" curr_x = x[i-1]\n",
	" prop_x = curr_x + np.random.normal(0, 1)\n",
	" a = target(prop_x) / target(curr_x)\n",
	" if np.random.uniform() < a:\n",
	" x[i] = prop_x\n",
	" else:\n",
	" x[i] = curr_x"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 7,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"0.5840424474631826"
	]
	},
	"execution_count": 7,
	"metadata": {},
	"output_type": "execute_result"
	}
	],
	"source": [
	"N = quad(target, -np.inf, np.inf)[0]\n",
	"N"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 8,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"[<matplotlib.lines.Line2D at 0x194f8310>]"
	]
	},
	"execution_count": 8,
	"metadata": {},
	"output_type": "execute_result"
	},
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYXHWV8PHvqeruJN3pTtJLNhLoAAEJiDBGFlEGFccwanAfcJlxRmV8X2HcxpkgyoyoMwz4uiC4oDI6iDKIiIFEA8ouAZNAgISAJDEknXSS3qv37qo67x/3VnWl00t19626S5/P8/CkqvrWvccyOXX6/JYrqooxxphoifkdgDHGGO9ZcjfGmAiy5G6MMRFkyd0YYyLIkrsxxkSQJXdjjIkgS+7GGBNBltyNMSaCLLkbY0wElfh14draWq2vr/fr8sYYE0pbtmxpVtW68Y7LK7mLyCrgW0Ac+KGqXjvs598A3uA+LQfmq+rcsc5ZX1/P5s2b87m8McYYl4i8nM9x4yZ3EYkDNwFvBhqATSKyVlWfzxyjqp/OOf4K4MwJR2yMMcYz+fTczwJ2qupuVR0AbgcuHuP4S4GfexGcMcaYycknuR8D7Mt53uC+dhQROQ5YBjww9dCMMcZMVj7JXUZ4bbR9gi8B7lTV1IgnErlMRDaLyOampqZ8YzTGGDNB+ST3BmBpzvMlwIFRjr2EMVoyqnqzqq5U1ZV1deMO9hpjjJmkfJL7JmC5iCwTkTKcBL52+EEicjIwD9jobYjGGGMmatzkrqpJ4HJgA7ADuENVt4vINSKyOufQS4Hb1W7tZIwxvstrnruqrgfWD3vt6mHP/927sIwxxkyFbT9gjDER5Nv2A8YUUv2addnHe659q4+RGOMPq9yNMSaCLLkbY0wEWXKPioYt8MBXoC2vPYWire1lPlNyB2fKS35HYoxvLLlHQfNL8JO3wyPXw0/fBYN9fkfkn4EeuPWd/FPJ3fy87CsskcN+R2SMLyy5h1D9mnXZ/wB49OuAwsXfgZadsPWnvsbnq623Qesu/nnwH0kT41Mld/kdkTG+sOQecqes+SV9W3/Brb2vhTM/AAteCU9P4+T+1P/AoldxZ+ovuTv1Wi6KPQkD3X5HZUzRWXIPufNjzzJTBlmfPtt54ZXvgQNPQ2K07X8irPMQHHwWTn0XAGvT51Eh/bDz9z4HZkzxWXIPufNjz5HQWfwx/QrnhRPf5Py560H/gvLLnx92/jz+AgC2pE+iW2fAnx/xLSRj/GLJPeT+IvYSW9MnkiLuvLDgNCivhZf/4G9gftj9MMyaBwtPB2CQEjalXwF7HvU5MGOKz5J7iFXQy0myj6d1OeAOtF65ngc6lzitmelm70Y47jyIDf21fiq9HJpehP5OHwMzpvgsuYfYq2K7iIs6CSzHc3o8NL0wvQYS+zuhdRcsOuOIl7dpPaBwcJsvYRnjF0vuIXa67Abg6fQJR7z+bHoZaBoan/UjLH9kkvei0494+bn0MudB4zNFDsgYf1lyD7HlsQYatZoEs494/dlMsp9OrZmDzzl/LnzlES83MQ9mL7DkbqYdS+4hdpI08FL66HuVNzEXZlU7rZnp4uAzzkBy5aKjfvRgx0Kef/oPR+wUaUzUWXIPKSHNiXKAP+mSkQ+oO9nZlmC6OPwCLFgBcvT93HfqMSyTRoS0D4EZ4w9L7iG1VJqYJQOjJ/fa5dD8YnGD8osqtLwENctH/PFuXcQsGWAxLUUOzBj/WHIPqeXSAMDOEdoyANSeDD0t0D0NElpPC/R1OF9oI9idXgzA8bHGYkZljK8suYdUvRwCYJcuHvmAupOdP5v/VKSIfNSy0/mz5sQRf7xbnT788WLJ3UwfltxDaqkcJqHldFAx8gGZKnY6JPfM2MIoyb2JOSS0nONlGu63Y6YtS+4hdawcZq/OB44eQASgaglIHNr3FjUuX7TshFgpzD12lAOE3bqQZXKwqGEZ46e8kruIrBKRF0Vkp4isGeWY94nI8yKyXUR+5m2YZrih5D6KeAnMOQbap8GdmVp3w7x6iMVHPaRB57NEmooXkzE+Gze5i0gcuAm4CFgBXCoiK4Ydsxy4EjhPVU8FPlWAWI1LSLNUmtindaMeU79mHY+3zGbLM1uLGJlPOvaNUbU7GrSOY6QZ0jYd0kwP+VTuZwE7VXW3qg4AtwMXDzvmY8BNqtoGoKp2b7MCmk87M2SQfWNV7sA+nc/S6VCttu/NI7nXMkOS0HWoSEEZ4698kvsxwL6c5w3ua7lOAk4SkT+IyBMismqkE4nIZSKyWUQ2NzVNg6RTIMe69wUdsy0D7NM65ks7DPYWIyx/DHQ7UyHnLh3zsAatdR507BvzOGOiIp/kPtKInQ57XgIsBy4ALgV+KCJzj3qT6s2qulJVV9bVjd5SMGPLN7k3ZNo27RFOaJn/bXOPG/Owoc9iGgwwG0N+yb0ByC2LlgDD55Q1AL9W1UFV/TPwIk6yNwVwbOwwaRX2j9FzB4Z68lEeVM0k6zljV+77M5V7lD8LY3Lkk9w3ActFZJmIlAGXAGuHHXM38AYAEanFadPs9jJQM+QYaeYQ8xikZMzjsj35tj2FD8ovHW5yH6fn3stMWrTSKnczbYyb3FU1CVwObAB2AHeo6nYRuUZEVruHbQBaROR54EHgc6o6Dda9+2MhLTRq9bjHNTGHfi2NdrXavhfiZc62vuNo0DpL7mbaGLv0c6nqemD9sNeuznmswGfc/0yBLZJWXtCx2xAASowDWs2yjv1FiMon7ftgzpIjbq03mgat5VVRHn8wJoetUA0bVRZKKwe1Jq/DD2oNJCK87L5977j99owDWuvMltHh8wGMiR5L7mHT10GF9OfVlgE4yDzojHByT+x3Kvc8HNR5kOyD3rYCB2WM/yy5h41bhR/MN7lrNSQao7kyM51yFiWNcPelkRzKfGadtseMiT5L7mHjJve8K3ethvSgs9AnarqbnBuBV+WX3A/qPOdBlH+TMcZlyT1sEs7g6IQq95z3RUqnuz97npX7QTKfhe3rbqIvr9kyJkASB0ircIh5eR0+VK02AmcULq4iytzo+sLYFn5YBlQuzOt9h4/4LIyJNqvcwyaxnybmkMzze7kxM6smgpX7AnEHRitHuRvVMAOUwqxqS+5mWrDkHjaJA3n32wGamePctCOCrYj50gYSg4oJ7FNUtTiSn4Uxw1lyD5vOg0PthTykiTltiwjOdV9IG1TMd25Mkq/KRTagaqYFS+5h03WIJj1qw82xRTShLZC2vGfKZFUtsqmQZlqw5B4mKWdKYxNzJva+qkWRbEUskLa8Z8pkVS6CrsPOZ2lMhNlsmTDpbgZ04pX77IXw50cLEpKfFkjrETNlMrNoxlS5CFBn8VOeK1uNCSOr3MPEvUVck06wcp+9APraIdlfgKD8UcYg1dKV90yZrMyXgd1uz0ScJfcw6XLuwDTxyn3+Ee+PgvnS7jzIc457VgQ/C2NGYsk9TLrd5M5Ek7u713mEEtoCWp0HE+y5v/bb2wD41//5vdchGRMoltzDZNJtmUy1Gp1WxNACpolV7i1UAVBHu9chGRMoltzDpOswzJhDP2UTe1+2co9Ocq+VDufBBJN7P2V0aPnQ+42JKEvuYdJ1CGZPYDVmRmYFZ4TaMjWSIK0Cs/Jf0JXRrHMsuZvIs6mQYdJ1OK97hR6lpMzZUyVClXsNCVqpZOXnfzvh9zYxlzpL7ibirHIPk65DQ/3ziZq9IFrJXRK0auWk3tusc6jFkruJNkvuYTLZyh2cL4UItWWqpZOWiQ4su5p0DnViA6om2vJK7iKySkReFJGdIrJmhJ9/WESaRGSr+99HvQ91mhvshf6EVe6uWjqyM18mqknnUCW9MNjncVTGBMe4PXcRiQM3AW8GGoBNIrJWVZ8fduj/qurlBYjRwFDVPdnKvXKBcw5VEPEuLp/USIKW9CTbMpm9eboPw9xjPYzKmODIp3I/C9ipqrtVdQC4Hbi4sGGZo0w1uc9eAMle6O/0LiaflJBkrnRPui3TnHlfhNpUxgyXT3I/BtiX87zBfW24d4vIsyJyp4gs9SQ6MyTTUplKWwYikdDm4XxBTb4t467wjcBnYcxo8knuI/0Or8Oe3wPUq+rpwO+An4x4IpHLRGSziGxuamqaWKTTXSa5V0w2uUdnlWqtJABo0ckl96HKPfyfhTGjySe5NwC5lfgS4Ig7P6hqi6pmthz8AfDqkU6kqjer6kpVXVlXN4nFONNZ12FAoKJ2cu+viE5yr55ics9W/N1WYJjoyie5bwKWi8gyESkDLgHW5h4gIrm7N60GdngXogG47febaNHZ1F913+ROkFml2t3sXVA+qcFN7pNsywxQSrtWWFvGRNq4s2VUNSkilwMbgDhwi6puF5FrgM2quhb4JxFZDSSBVuDDBYx5WqqWxFA7YTLKqwGBnvAn96m2ZcBpzcyNwG8xxowmr+0HVHU9sH7Ya1fnPL4SuNLb0EyuaumkdQrJjFgcymsi0YqolgSDGidB+aTP0UIVJ/a0eBiVMcFiK1RDooYELUxuXndWRW0kkntmXxmdwl/fVq2MRIvKmNFYcg+JaklMrXIHp+/eHf5qtdaDz6JVqyLRojJmNJbcwyCdYi7dtE61co9QW6Z5ism9hUroaYV0yqOojAkWS+5h0NNKTHRKA4iAU7lHoFp1WlQeVO4o9LZ5E5QxAWPJPQzchDz1tkytk8xSgx4E5Z8ar9oyYH13E1mW3MPATUBTbstkFkCFeZbIYB+V0jvltkxzpvKPwG8yxozEknsYuAnIk7YMhLtazfwW40lbhnB/FsaMwZJ7GHR71JYpdyv3MA+qdnvzRdeSuYuTVe4moiy5h4HbRmlj9tTOk6ncw9yW8Si5t2VaXBGYGmrMSCy5h0F3M206mxTxqZ2nIgqVuxP7VGfLJCmBmXOscjeRZck9DHqaJ30z6CPMnAsSD3dy92r8AZw2lfXcTURZcg+D7uYpV6oAxGLuQqYQJ7TuJvq1lC5mTf1cFbVWuZvIsuQeBj0ttHlRuYO7BUGIE1p3s7vHjgf3gS2vtZ67iay8doU0PutupkVPnfTb69esyz7ec0rIq9XuZm9aMgAVNbB/szfnMiZgrHIPunQaelqmPK87K+w7Q3Y3TfrG2Ecpr3VmDunwu0YaE35WuQdY/Zp1zKGLZ2ampj7HPSMSbZljvTlXRS2kk9DXDrPmeXNOYwLCKveAq8nedcijnnt5LfQnINk//rFB1NPsbeUO1nc3kWTJPeCq3fuFetqWgXBW7wPdMNjjbc8dwj0GYcwoLLkHXI10Ah5sPZCR3TwshAnNHSuY8gZqGeUh/qIzZhyW3AOu2uu2THbzsBAOqrrtkyndKDxXmL/ojBmHJfeAq8at3D1ry2SSewj7zJmtB7z6LcYqdxNheSV3EVklIi+KyE4RWTPGce8RERWRld6FOL3VSIJOncUApd6csNztM4excvdou9+s0plQNjvcG6kZM4pxk7uIxIGbgIuAFcClIrJihOMqgX8CnvQ6yOnMuTG2Ry0ZcDbLipWGM7lnK3cPP4+wb8dgzCjyqdzPAnaq6m5VHQBuBy4e4bgvA9cBfR7GN+1V0+ldpQogEt49VbqbobSCXmZ6d86wfhbGjCOf5H4MsC/neYP7WpaInAksVdV7PYzN4LRlPK1UwV2lGsKE1t08NH3RK7YzpImofJL7SDs0Zddri0gM+Abw2XFPJHKZiGwWkc1NTSFsC/igWjq9mwaZUR7SLQi6m4YGhL1SUWs9dxNJ+ST3BmBpzvMlwIGc55XAacBDIrIHOAdYO9KgqqrerKorVXVlXZ3H/0gjSakm4W1bBsK7BUF309AMF69keu62v4yJmHyS+yZguYgsE5Ey4BJgbeaHqtqhqrWqWq+q9cATwGpVte32pmg2vcyQZGHaMmGsVntaClO5p/phoMvb8xrjs3GTu6omgcuBDcAO4A5V3S4i14jI6kIHOJ1Ve706NaO8xklmg73enreQVN22jNeVu811N9GU166QqroeWD/statHOfaCqYdlAGrcfWVavFpun5FdyNQMc5eOfWxQ9CcgNeB9cs+uUm2B6mXentsYH9kK1QDLbD3geeUexhtlZyprr9syVrmbiLLkHmDZtkwhBlQhXH33bHL3ui1T7fxpc91NxNjNOgIss6+M5wOq2S0IQpTQMr9llNcC+z05Zf2adVTQy/aZhOuzMCYPVrkH2DzppE9L6WWGtycO486QPYVpy3Qzk34tDddvMcbkwZJ7gNWQoIUqRl5HNgUzKiFeFq5WROaLyOu2DOIMWFtyNxFjbZkAmyedtHnckqlfsw6AjTNmsyhMrYjuFphRBSUe/xaDM2C9OEyfhTF5sMo9wGqk09sdIXO0aFW4+syFmOPuatXKcP0WY0weLLkHWHW2LeO90CW0Quwr42qlMlxfdMbkwZJ7gBWiLZPRzJyQDai2eL+vjKtVq6CntSDnNsYvltyDKtlPlfR6d0u5YVq1Mly32itgW6ZFq2CgE5L9BTm/MX6w5B5UbiXZ5vXWA65WrYLBbhjoKcj5PZVOu3u5F7AtA9aaMZFiyT2o3H645wuYXM2ZXn4Y+u597aCpwg6oQjg+C2PyZMk9qNx514XquWf3qwlD371Q+8q4hj4LS+4mOiy5B5WbaAo1Wybbyw9D3z279YDHt9hzZdsyNqhqIsSSe1Bleu6FmuceprZMgbYeyMh+0YXhszAmT5bcg6qnmbQK7cwuyOmzfeZQtGUyWw8UJrl3UAESs7aMiRRL7kHV00I7FaQL9H9RF7MgPiMcCS0TY2Z7Xo8pMZhVbZW7iRTbWyaoupu9v0nHEcSZfRKS5N6msznzqvsKd42wfBbG5Mkq96DqaRka6CuUitpwVKvdTQVbzJVVXmsDqiZSLLkHVU9LgSt3nIQWhmq1u7lgs4ayKmrC8UVnTJ4suQdVTwutWpjB1KywtCJ6motTuYfhszAmT5bcg0jVbcsUulqtC0e1WpS2TA30tkE6VdjrGFMkeSV3EVklIi+KyE4RWTPCzz8uIs+JyFYReUxEVngf6jTS1wHpZMHmuGeV18BgDwx0F/Y6U5FOQU9rEb7oagG1vruJjHGTu4jEgZuAi4AVwKUjJO+fqeorVfUM4Drg655HOp24Ww8UvFrN3ks1wNV7TyugNBejcge73Z6JjHwq97OAnaq6W1UHgNuBi3MPUNVEztMKQL0LcRpyk22hdoTMymzEFeTWjLuAqeCDy2H4LIyZgHzmuR8D7Mt53gCcPfwgEfkE8BmgDHijJ9FNV9nKvdBtGTehBblyd5N7s84p7HXC8FkYMwH5VO4ywmtHVeaqepOqngD8K/CFEU8kcpmIbBaRzU1NIVj27he3eixatRrkhJZJ7oXuuWfbMgH+LIyZgHySewOwNOf5EuDAGMffDrxjpB+o6s2qulJVV9bVFWafkEhwK/eiLGKCYO8vU7TK3U3uYdgl05g85JPcNwHLRWSZiJQBlwBrcw8QkeU5T98KvORdiNNQdzOUltPHjMJep2w2lMwMdrXa3QSxEhKUF/Y6JWUwY44NqJrIGLfnrqpJEbkc2ADEgVtUdbuIXANsVtW1wOUiciEwCLQBf1fIoKOsfs06vlb6LOfECpzMAETcxTsBTmhdh6G8Fu0pwpIMW6VqIiSvjcNUdT2wfthrV+c8/qTHcU1r1SSGtuQttIqagLdlmmF2HRQj59oqVRMhtitkAFVLZ+EXMGUEeJVq/Zp1/KrsJTp1VnEuWF4DHfvGP86YELDtBwKomkThN8rKCHi1WksHzRR4MDWjoibQn4UxE2HJPYDmSVcRK3c3uWsw153VSKLwM2Vwfkv47qYOBjqbAvtZGDMRltwDpoxBKqW38AuYMipqIdkbyP1lZtFHufQXfhsGV4tWUSYp6E+Mf7AxAWfJPWDm0QkUYeuBjMz+MgHsu9dKB0DRWlRD95UN3mdhzERZcg+YGnGqxoKvTs3ILrsP3nTIWpzPohhtGchZNGZz3U0EWHIPmHniVO7FmwoZ3FWqmS+6gu8I6cp+oVrlbiLAknvA1LhtmYJvPZAR4N0Qs22ZolXubnIP4GdhzERZcg+Y6mxbpvDJvX7NOlb81xYArr3z0YJfb6Jq3LZMsXru2UFsa8uYCLDkHjA1kiCpMdop8P1TXT3MoFfLqHbbQUFSKx0ktJwBSotyvV5m0qtl1pYxkWDJPWBq6aCFKrRo/9cILVRl+9tBUisdReu3Z7RSaZW7iQRL7gFTK4mi9ZgzWrUy2wIJkhoSxVud6mrVSqvcTSRYcg8YP6rVJp2bHbwMkhpJFG0BU0arVtmAqokES+4BU1PMvVRcTTqHOmkv6jXzUSsdRU/uLVRZW8ZEgiX3gPGjWm1iLrV0QDpV1OuOKZVkHl0+tWUsuZvws+QeJAPdVBRxL5WMwzqXuGiwKtaeFmKiRVudmtGqlTDYDYO9Rb2uMV6z5B4k2ZtBF7stM9d50HWoqNcdk/tZFL3njq1SNdFgyT1I3IRS/AFV98uk05J7dvGYDaqakLPkHiRdh4HiLbfPOEwAK3c3lmxsRZL9MglSi8qYSbDkHiSZtkyRk3v2ekFK7p0HATis84p6WWvLmKiw5B4k3U7lXrRNw1y9zHTuUxqk5N51iC6dSQ8zi3rZ7BddAHfJNGYi8kruIrJKRF4UkZ0ismaEn39GRJ4XkWdF5Pcicpz3oU4D3c0ktJx+yop+6cM6N1jJvfPg0FhAMS/LLCiZGazPwphJGDe5i0gcuAm4CFgBXCoiK4Yd9jSwUlVPB+4ErvM60Gmhu6nog6kZzczJ9vwDoesQhyluS8YhMHt+sD4LYyYhn8r9LGCnqu5W1QHgduDi3ANU9UFV7XGfPgEs8TbMaaLrcNG2tx2uKZCVe3EHU7NmLwjWZ2HMJOST3I8B9uU8b3BfG81HgN9MJahpq7u56IOpGYd1brCmQnYdcmLyw+wFVrmb0MsnucsIr+mIB4p8EFgJXD/Kzy8Tkc0isrmpyQasjtLdVPR53RlNOhcGOmGg25frH6G/Cwa6/EvuFXVWuZvQyye5NwBLc54vAQ4MP0hELgSuAlarav9IJ1LVm1V1paqurKurm0y80ZVKQk8LLUVenZrRlLluECrWzBx3Pyv3nlZIDfpzfWM8kE9y3wQsF5FlIlIGXAKszT1ARM4Evo+T2AOQHUKotxVQX2aIQE4iDULFmpnj7suAKs6AKmpz3U2ojZvcVTUJXA5sAHYAd6jqdhG5RkRWu4ddD8wGfiEiW0Vk7SinM6PJrk71sS0DwUjuXZkFTD5W7hCMz8KYSSrJ5yBVXQ+sH/ba1TmPL/Q4runHXcDk14DqUHIPwC9enQFoy0AwPgtjJslWqAZFthXhT0JrpRIkno3DV10HIV5WtJuEH2X2fDcOq9xNeFlyD4pOf1sRaWJOUgtCcu885FbPI03UKgJL7iYCLLkHRdchmFFFb5H3UjlC1WLoPGoiVPF1HRxqjfihdBbMCNiKXWMmyJJ7UHT6nNDASe6JACT3zkNQudDfGGbPt8rdhJol96DoPOh/Qqs6JhjJPbHficVPsxfYzpAm1Cy5B4XfrQhwKvf+BPQl/IuhL+HEULXYtxDq16zj3t1Jdv15N/Vr1vkWhzFTYck9CFSD0YrIVMudjf7FkPnNYY6/e8816RzqpN3XGIyZCkvuAfDKK++EZC9feaTN30Ay1XJiv38xJBrcWPxtyzTpXKqkl5mMuJOGMYFnyT0A5ouT1H1btON63fdeBOBzt/i4qWemcvexLQNDe+3USoevcRgzWZbcA2C+++t/k08LmDIy9ytdSKt/QXTsBwQqF/kXA0MrdudjrRkTTpbcA2A+wajcByilSatYJP4l99t//wSHdQ71X7jftxgAGrUagIU+fhbGTIUl9wDIVO6H1KddEHMc0mpfE9piaeGAm1j9lEnufn7RGTMVltwDYL6006Mz6GKW36HQqNUsFP8GdhdJK41a49v1MxJU0KtlVrmb0LLkHgALpM1tyfi0l0qOg1rNQmnx7foLpZWDAajcQWjUaqvcTWhZcg+A+dLu226QwzVqNdXSBYO9xb94XweV0suBAFTu4HzRLbDkbkIqr/3co2y0FYh7rn1r0WJYRAvP6AlFu95YDmit86CjAWqXF/fi7jTIYFTu0Eg1Z8sLfodhzKRY5e63dJqF0hqYanWfuve2bdtT/It3OIungvJZHNJqFtAG6bTfoRgzYZbc/dbTzAxJBmIQEWCfunuZ+5Lc9wI5vz34rFGrKZUU9Ni9VE34WHL3W4ez3L4xIK2IJubQr6XQ/nLxL962h34t4aBfN8YeJtse8nM7BmMmadr33H2XCFYrQomxT+s4sc2P5P4yDVqHBqTmyH7hJg7A4jOLcs3cMaBijvuY6AnGv6LpzB1EDEpbBty+ux9tmbY9Q22hADiUm9yNCRmr3Ecwnzb4w7egfR/MfwW88n0ws6owF+tooF9LaaFA55+EfTof2v9Y/Au37WGfriz+dUfRTBWDGqe0kMm9vwue+Tk0vQBzj6OWWprdTcuMmYq8kruIrAK+BcSBH6rqtcN+fj7wTeB04BJVvdPrQIvl4thjfLX0Fri/z7mPZn8HPPoNuOQ2WHyG9xdM7Hd//fd/AVPGPq2Dvg7obYdZRZp/39sOfe3sDVDlrsQ4xDyWFCq5H9gKP7/E2T/f/bv28IwZXDn4UdamzyvMNc20MW5bRkTiwE3ARcAK4FIRWTHssL3Ah4GfeR1gMV0Sf4BvlX2H7VoPl2+BNS/DP2wAEbj1HdC80/uLduwPVEsGGEqwxRxUda8VpLYMuH13d9DbSxdc+UMS37+IhsQg/MN9cOVeuHwL23QZN5TdxHvjD3l+TTO95FO5nwXsVNXdACJyO3Ax8HzmAFXd4/4stBOCz41t56slP+LB1Kv4+OCn6f/ai4Czv/lS+Qy/LvsiB254N+8cuIZB92PzZMArsZ8DHDf183ioIXc65KJXFeeibcFM7g1ax2va93p70tQg3y39FkliXDLwRR479mzn9doT+eDA5/lR6fVcW/ID2P1WOP4Cb69tpo18BlSPAfblPG9wX5swEblMRDaLyOampuDcfLiGDr5d+m1262I+MfhagDdEAAASkklEQVRJ+ik74uf7dAFXDn6M02J7+Nv4fd5dOJWExIHAVe5DC5mKWLm7A7hBasuA+1kkGiA16N1Jn/w+p8T2cuXgx2jIfNauQUr4+OCn2a2L4a7LoNvm2JvJySe5j9QM1slcTFVvVtWVqrqyrq5u/DcUyRdLb6WSHv7v4CfpYeaIx2xIv4ZHUq/kipJfMYcuby6caABNBS6hJaiAmXOKNmOmfs06bvvNQ7TpbDopL8o187VP54Omef0XbvXmZtndLfDQf/JA6gw2pEcePO5hJlcMXuGMQ6z/3NSvaaalfJJ7A7A05/kSIDJzw86PPcM74o/z3dRqXtKxb8r81eQHqKKHj5V48I8csskzaK0IAGpOhJYCjDGM4nhpZJf6e2u9kWQq6yXiUQX95HdhoIv/SL6fsQbRX9Bj4fWfhe13wa4HvLm2mVbySe6bgOUiskxEyoBLgLWFDatIkv18ueS/2ZVexHeTq8c9/EU9lvvSK/lg/HfMom/q18+0ItIBTO61J0HzS0W73AmxA+xKBy+5Z1pUS+XwlM912ppfkHj4Jn6Teg07xykkAHjdp6D6BFj3WUgOTPn6ZnoZN7mrahK4HNgA7ADuUNXtInKNiKwGEJHXiEgD8F7g+yKyvZBBe2bzLRwXO8yXkn97VJ99NDcn38pc6ea98Yenfv22PRAroZFg9dwBZ0fIzgPQ31nwS1XRzXxpZ5f6e9/UkTRqDUmNsVSmPkb07vijVElPXoUEACUz4KLroHU3bL5lytc300te89xVdT2wfthrV+c83oTTrgmPvg54+DoeS53KI+nT837bU3oSW9PH88H470DVmSY5WW17YO6xpHsCuFC49iTnz+Y/wTGvLuiljpdGAGcQMWBSxGnUmqlX7qpcGn+AZ9LH82ye2zs7PX7lp6Wn8rqH/wvOuNQZCzEmDwHMKoVXv2YdN3zln6C3lWuTlzLRBUT/m3oDJ8X2c/Hnb6B+zbrJD7S17YF59ZN7b6HVnuz8WYTWzAniDOEEsecOTt99ypV7wyZeEdvHz1JvmuAbxfk72tvqrJo2Jk/TMrnX0cZH479hbepctunxE37/Palz6dUy/maqC02CnNyrl0GsxKncC+z42AEGND40BTNg9nmR3Lf8mC6dyT2pcyf81m16PJz2Htj4HUg0Ti0OM21Mi71lhlfWXym5i1KSfC35vkmdr4ty1qfP5u3xjVyT/BB9zJj4SXrbobctsMm9/qr7+H1ZHS899AirJlpsTtAJ0sjLupBkQP867tGFzJeHqWCStx4c7IXnf829qXNGnGqb129+b/wCPP9rePhaeLtV8GZ8065yP04O8jfxh/h56o3s1QWTPs+dqfOplF7eFHt6cifIzCEPaHIHp02yXLxfej/cybKXl3RS6+KKIjPQmxkbmLCX7oOBLtamXzv5IKqXwcq/h6duLcw2GCZypl1y/3TJnSSJ8+3kO6Z0nifTp3BY5/K2+MbJnSDT7sj0tgPoeT2OZXIQBroLd5G+BMtih9ieri/cNaYoMxaQGRuYsG13QUUdT6ZPmVog538OSmbCA1+e2nnMtDCtkvsrZC+rYxv579QqmqZ4t580MdalzuYNsa2T+nX9xjvuJakxln/9xSnFUUjb0suIicLBbYW7yCHn3M9rsPbXybVXF5DUGCfEJpHc+7vgTxtgxcWkiE86hvo166j/yia+1fsWeP5u2P/UpM9lpofIJvfMLJbcfuZnS35BF7P4XvJtnlzj3tQ5zJRBLoxtmfB7T5QD7NGF2U3IgihbTR98tnAXaXz2yGsF0ACl7NX5HD+Zyv1Pv4VkL5z6Lk9i+UHyr6G8Bn7/JU/OZ6Irssl9uL+QP/Hm+Ba+l3wbCWZ7cs6ndDkHtJq3xZ+Y8HuXSwM7A9xnBmikmhathMathbvIwWdp0ioOU6R94ydply7mhMn03LfdBZWL4NiJz5IZSRfl8Pp/ht0Pwa4HPTmniaZpktyVfy65gyat4sepVR6eNcb61Nn8ZewZZ/ZLvpIDHCeHAj2I6BCnom4sbOX+fLqeIN2sZCS7dDHLpBHSqfzf1NcBO++HU98JMQ//qb3mIzDnWPjdv0M6tLtsmwKbFsn9jbGneW38eW5MvnPUXR8n697UuZRJCl6YwEKm1l2USJqdAdxLZbjtWg+Hd0Cy3/uTD/ZC0w7nGgG3WxcxQ5IT2ynzhfWQGvCsJZNVMgPe8Hlo3MonvvBvU1tIZyIrUsl9pD57GYN8seRWdqYXc9uEVweOb6uewL50HWz7Zf5vanIGUfPaPMpnT6WXQ3oQ9k98XGFc+7dAOsnm9Enen9tjL6SPdR4cmsDg8va7nAp7SQHuC3v6+2D+Cj5bcgclJL0/vwm9SCX3kXw4/luWxQ5xTfJDBVokI9yTPtfpgXa35PeWg8+R1Bg7A7rcPtcf068ABPb8wfuT73WmkW4JQXJ/UZeS1Fj+LaqeVmer3lPfMbX9h0YTi8Obrub42EE+FL8fGLm4MdNXpJP7fNq4ouRufpc6k0fShbtd3D2pc0FTsOPX+b2hcSsv6ZK8d6L0UwezYcFp8PJj3p/85Y1Qd4pzjYDrp8wZAM935tCOtZBOwmnvLlxQJ63iwdSr+GzJL1hEnoWFmTaCOw9vypSvlt5CKUm+kvxgQa+0Q49lZ3oxTb/+AZfeuWDse6uqwoGneS59WkFj8lT9ebDlJ86e4iUefSElB2Dfk3D63xx5E8cA267H8Yp8K/fn7oSaE6m/oQHY72kcuZX5Evl77i/7F/699Cf84+BnPL2OCbfIVu6rY4/z5vgWvpZ8H3sKvk+4cE/qXM6O7WA+bWMf2tEAPS08O4kNy/zyj4/OgmQv7/nit7076d6NMNAFJ17o3TkLbHt6GXQdHHfzrrPW/JT0nx/jm4dOp9CzgBp0Pt9Mvpu3xDezKvbHgl7LhEskk/siWvhS6U94On0it6QuKso170mfS0yUt443533/ZsBZ/RkWj6dPZUDjXBj3cFXkS/dBvAyWne/dOQtsS3q582Dv2FtOvC3+BDFR1qamsJfMBPwodRHPppdxbekPWIzdUNs4IpfcS0lyY9kNlJLkM4P/h3SR/ifu1sVsS9ezery9ZvY8BqUVbAvB9L+MTsp5PH0aq2KbnLbSVKnCi+vhuPNgRvD77RnbtR5Ky8dN7qvjj7MtXV+0m48kKeGKwSuIk+aGshtt9owBIpfclS+V/JhXx17iXwcv489Fvm3b2tS5nBnbOeoNLurXrOOlJ3/DQ30nBnZ729FsSK+kPnYIGp+Z+skaNju3jivkYGMBJCnh0b7j2fHEb0efkXJ4B2fEdnF36ryixvayLuSqwY+wMvYnrin5sTdfwibUIpXcPxm/i/eXPMB3kqtZlz6n6Ne/K3U+gxqHLT8e8ed1tLM8tp/H0yuKG5gH1qXOpk9L4amfTP1kz/yMPi3ltDtmhm7a3pPpUzhZ9lHHKCuSt/yEAY1zV+r1xQ0MWJt+Ld9Jrub9JQ/AI9cX/fomWEKf3J15vffy7S/8HZ8u/SW/SJ7Pdcm/8SWWZuawIb0Stt7mrL4c5k1uz/qx9CuLHdqUJZjNPalz6dr0M05fc8fkk3J3MzxzO/emz3X2SQmZ36VfTUyUC+MjLOoa7IVnfs6G9Gtopar4wQHXJ9/HL1Ovgwe/Cg981Sr4aSz0yb2CXr5ZehNXlNzNz5NvYE3yY/i5T8ltqQudOyw9c/tRP3t7bCO70osCvb3tWG5JXUQ5/Xy85J7Jn2TjjTDYy3eTb/cusCJ6QZfycno+b4ltPvqHW38Gfe38NPnm4gfmUmJ8bvDj/G/yAnjkOu6++q+dbYfNtJNXcheRVSLyoojsFJE1I/x8hoj8r/vzJ0Wk3utAR/TS/dxb9nneHtvI1wbfy5XJj05pz2wvbEyvgGNeDY993ZnLndGxn3Niz3Nv+hyCvknWaHbocdydPo+PxH/DSTLxyel/deX3GHj0Bu5KnceuwG+aNhphXfocXhd7Dtr3Dr2cHIDHvglLzuJJfYV/4eHca2BN8qN8bfC9vD22Eb7/etj5O19jMsU3bnIXkThwE3ARsAK4VESGN40/ArSp6onAN4D/8jrQrOSAs0DkR2+B296DIrx/4AvcmHonwUiaAhdc6fzDf/yGoZc33ogi/CL1l/6F5oH/HHw/Ccr5Xuk3IDGB/c0Tjfyg9P+RoIIvDxZ2UVmh3Zp8M4rAxpuGXvzDt6BjL1ywhiD8PVRi3Jh6Jx8YvAokBj99t/NvZvuvRmwZmujJZ8rGWcBOVd0NICK3AxcDz+ccczHw7+7jO4EbRURUC9Dwe+Q6Z7BoXj2s+i/ecveiwN3wov5HA9xUehYX/v4/mbH4TGc+9x9v5pep82nQ+X6HNyVNzOXjA5/ix2XXwffPh/P/BZa/GeYe6+x3kpFKQqofOg86VePD11ErnXxg4CrafOpHe6WRGn6Veh3v2/RDOGU1DPbS/8C13J8+m8t/2Od3eEd4Ir0C/s/jsPm/4YnvwC8+TK+W8WT6FC644M2wYAVUHcNffnc7bTqbAUp54T8uPvL/y6DJTSs5j5d93hkHUoQ9//nWwuzpMwW541RjrmL3iIyXf0XkPcAqVf2o+/xDwNmqennOMdvcYxrc57vcY0ZdUbFy5UrdvHmEvuV42vfy4ev+h4fTr0IDPGQwl07uKLuGk2LO0vNd6UW8a+BLodhHJR8ny17+o/RHvDqWM+0zVgqxEiep65H7jG9On8SawY+GYifMfMyhi1+VORt3AbyYXsL7Bq4O9P+/MdK8NradN8We4rWx7Zwcb3T2RBpBUmMkiTOzdKKFk9I36JxTcl4DmFESyybjwdTQ34/S2NBxqbQe8d6YeFMfplVQIB6TobPLWI/dKEZ8PDld/UPrD2ZffD38xd9O6jwiskVVx91qNJ/k/l7gLcOS+1mqekXOMdvdY3KT+1mq2jLsXJcBl7lPTwaG30C0FgK5xM7imhiLa2IsrokLamzFiOs4Va0b76B8vpYbgKU5z5cAw5utmWMaRKQEmAO0Dj+Rqt4M3DzahURkcz7fSMVmcU2MxTUxFtfEBTW2IMWVT19jE7BcRJaJSBlwCbB22DFrgb9zH78HeKAg/XZjjDF5GbdyV9WkiFwObADiwC2qul1ErgE2q+pa4EfArSKyE6div6SQQRtjjBlbXqMlqroeWD/statzHvcB7/UgnlFbNj6zuCbG4poYi2vighpbYOIad0DVGGNM+AR3LqExxphJC1RyF5Evi8izIrJVRO4TkUDcQVpErheRF9zYfiUic/2OKUNE3isi20UkLSK+jtKPt02FX0TkFhE57K7HCAwRWSoiD4rIDvf/w0/6HROAiMwUkT+KyDNuXF/yO6ZcIhIXkadF5F6/Y8kQkT0i8pybuyaxgMd7gUruwPWqerqqngHcC1w93huK5H7gNFU9HfgTcKXP8eTaBrwLeMTPIPLcpsIvPwZW+R3ECJLAZ1X1FOAc4BMB+cz6gTeq6quAM4BVIlL8PbRH90lgh99BjOANqnpGmKZCFo2qJnKeVpBZ2uYzVb1PVTPLy57AmesfCKq6Q1WHLwbzQ3abClUdADLbVPhOVR9hhHUXflPVRlV9yn3ciZOwfN9RTR2ZrSRL3f8C8W9RRJYAbwV+6HcsQReo5A4gIl8VkX3ABwhO5Z7rH4Df+B1EAB0D5G4V2UAAElVYuDupngk86W8kDrf1sRU4DNyvqoGIC/gm8C9AerwDi0yB+0Rki7sS33dFT+4i8jsR2TbCfxcDqOpVqroUuA24fOyzFS8u95ircH6Vvq1YceUbWwCMtPFGIKq9oBOR2cAvgU8N++3VN6qactujS4CzROQ0v2MSkbcBh1V1hDul+O48Vf0LnLbkJ0TE9zu/F307RVW9MM9DfwasA/6tgOFkjReXiPwd8DbgTcVefTuBz8xP+WxTYYYRkVKcxH6bqt7ldzzDqWq7iDyEM2bh94D0ecBqEflrYCZQJSI/VVXf95BW1QPun4dF5Fc4bUpfx8EC1ZYRkeU5T1cDL/gVSy4RWQX8K7BaVXv8jieg8tmmwuQQEcFZ3b1DVb/udzwZIlKXmREmIrOACwnAv0VVvVJVl6hqPc7frweCkNhFpEJEKjOPgb/C/y/CYCV34Fq33fAszgcUiKlhwI1AJXC/O9Xpe34HlCEi7xSRBuBcYJ2IbPAjDnfAObNNxQ7gDlXd7kcsw4nIz4GNwMki0iAiH/E7Jtd5wIeAN7p/r7a6VanfFgEPuv8ON+H03AMz7TCAFgCPicgzwB+Bdar6W59jshWqxhgTRUGr3I0xxnjAkrsxxkSQJXdjjIkgS+7GGBNBltyNMSaCLLkbY0wEWXI3xpgIsuRujDER9P8BtSd1JqHhvXsAAAAASUVORK5CYII=\n",
	"text/plain": [
	"<Figure size 432x288 with 1 Axes>"
	]
	},
	"metadata": {
	"needs_background": "light"
	},
	"output_type": "display_data"
	}
	],
	"source": [
	"plt.hist(x, bins=100, density=True)\n",
	"f = np.vectorize(target)(np.linspace(np.min(x), np.max(x), 1000)) / N\n",
	"plt.plot(np.linspace(np.min(x), np.max(x), 1000), f)"
	]
	}
	],
	"metadata": {
	"kernelspec": {
	"display_name": "Python 3",
	"language": "python",
	"name": "python3"
	},
	"language_info": {
	"codemirror_mode": {
	"name": "ipython",
	"version": 3
	},
	"file_extension": ".py",
	"mimetype": "text/x-python",
	"name": "python",
	"nbconvert_exporter": "python",
	"pygments_lexer": "ipython3",
	"version": "3.7.3"
	}
	},
	"nbformat": 4,
	"nbformat_minor": 2
	}