{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import quad\n",
    "from scipy.stats import norm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pdftarget(x, normed=1):\n",
    "    return np.exp(0.4*(x - 0.4)*(x - 0.4) - 0.08*x*x*x*x)/normed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(7.8521781788662155, 1.4333237621618567e-08)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f7328d1a150>]"
      ]
     },
     "execution_count": 3,
     "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 0x7f732ac08350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# integrate \n",
    "I = quad(pdftarget, -100, 100)\n",
    "print I\n",
    "\n",
    "xmin, xmax = -5, 5\n",
    "N = 100\n",
    "x = np.linspace(xmin, xmax, N)\n",
    "y = pdftarget(x, normed=I[0])\n",
    "\n",
    "plt.plot(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sampling(mu, sigma, nsamples=1000, xmin=-5, xmax=5):\n",
    "    I = quad(pdftarget, -100, 100) # integrate\n",
    "    normalize = I[0]\n",
    "    \n",
    "    # find good propos dist\n",
    "    x = np.linspace(xmin, xmax, 1000)\n",
    "    target_dist = pdftarget(x, normed=normalize)\n",
    "    propos_dist = norm.pdf(x, mu, sigma)\n",
    "    \n",
    "    xhighest_idx = np.argmax(target_dist)\n",
    "    yhighest = np.amax(target_dist)\n",
    "    c = 1.05*yhighest/propos_dist[xhighest_idx]\n",
    "    \n",
    "    # start sampling\n",
    "    xsample = np.random.normal(mu, sigma, nsamples) \n",
    "    target_dist = pdftarget(xsample, normed=normalize)\n",
    "    ysample = c*norm.pdf(xsample, mu, sigma)*np.random.random(nsamples)\n",
    "    \n",
    "    accept = ysample <= target_dist\n",
    "    \n",
    "    return xsample, ysample, accept, c, normalize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VEXbx/HvnUIgCZ0QeiehBEITQVoA6VIERIogCAIq\n6AtYqBKaiNhQihRFBQFRUHoJJfReAoEUeu8ECKGlzPtHgk/EkGyS3T2bzXyui+vJ7s6c+SUP3pzM\nOWdGlFJomqZp9sXB6ACapmma+enirmmaZod0cdc0TbNDurhrmqbZIV3cNU3T7JAu7pqmaXbIpOIu\nIs1FJFREwkXkk2TavSAi0SLSPtF7Z0UkSEQOichec4TWNE3TkueUUgMRcQCmAo2By8A+EVmmlApN\not3nwLpnDhEH+CmlIswTWdM0TUuJKWfuNYETSqlzSqloYBHQNol2A4E/gevPvC8mjqNpmqaZiSlF\ntzBwIdHriwnv/UNECgHtlFIziC/miSkgQET2icjb6QmraZqmmSbFaRkTfQsknotPXODrKKWuiIgH\n8UU+RCm13UzjapqmaUkwpbhfAoolel0k4b3EagCLRESAfEALEYlWSi1XSl0BUErdEJG/iJ/m+U9x\nFxG9yI2maVoqKaWenS0BTJuW2QeUEZHiIpIF6Awsf+bgpRL+lCR+3v1dpdRyEXEVEXcAEXEDmgLB\nyYQ09M/o0aMNz2Arf/TPQv8s9M/C9n8WyUnxzF0pFSsiA4D1xP9j8KNSKkRE+sV/rGY92yXR157A\nXwln5U7Ab0qp9SmNqWmapqWPSXPuSqm1gPcz7818Ttu3En19BqiSnoCapmla6ulbFBPx8/MzOoLN\n0D+L/9E/i//RP4v/sfWfhaQ0b2MtIqJsJYumaVpGICKodFxQ1TRN0zIYXdw1TdPskLkeYtK0dHsY\n/ZDjN44Tfiucx7GPiYmLISYuBg9XD8rmLUvp3KVxy+JmdExNyxB0cdcMExsXy5ZzW1h8bDGBZwM5\nd/ccZfOUpVy+crg6u+Lk4ISjOHIt6honbp/gdMRpiuQowitlX6GNdxvqFquLs6Oz0d+GptkkfUFV\ns7rzd8/z9a6vWRS8iMI5CtOpQidalG1BuXzlyOKY5bn9YuNiCb4ezPKw5SwPX87piNO8Xe1tPnjx\nAwpmL2jF70DTbENyF1R1cdes5uyds3y+/XMWH1tMn2p96Fu9L2XylEnz8c5EnOGb3d8w/8h82pdv\nz/B6wymVu5QZE2uabdN3y2iGehTziBEbR1B9VnXyZMtD+MBwvmjyRboKO0DJ3CX5rsV3hA8Mp0iO\nItScXZNxW8bxOOaxmZJrWsalz9w1i9p+fjt9lvfBJ78PU1tOpYB7AYuNdf7ueQauGUjYzTBmtJpB\nw5INLTaWptkCPS2jWV10bDTDNg5jwdEFfN/iezpU6GC1sZeFLuO91e/Rq0ov/P38cXRwtNrYmmZN\nurhrVnU96jqd/uhENuds/Nb+N/Jky2P1DNfuX6PLki44iAMLOiwgv1t+q2fQNEvTc+6a1ey/vJ8X\nZr9AnaJ1WNllpSGFHcDT3ZOA7gHUKlKL6rOqs//yfkNyaJpR9Jm7ZjarwlfRc1lPZr4yk/bl2xsd\n5x9/h/5N3xV9+a39bzQp3cToOJpmNnpaRrO4hUcXMmjdIJZ1XsaLRV40Os5/bDu3jY5/dGRK8yl0\n9ulsdBxNM4vkirt+QlVLt5n7ZzJu6zg29NiAT34fo+MkqV7xemzovoEWv7Xg9sPbvPvCu0ZH0jSL\nMmnOXUSai0ioiISLyCfJtHtBRKJFpH1q+2oZ09S9U5m0YxJbem6x2cL+VCXPSmzrtY1JOybx48Ef\njY6jaRaV4pm7iDgAU4HGwGVgn4gsU0qFJtHuc2BdavtqGdO8oHl8seMLtvbaSolcJYyOY5KSuUsS\n0D2Ahr80xNXZlS6VuhgdSdMswpQz95rACaXUOaVUNLAIaJtEu4HEb459PQ19tQxmedhyPt7wMeve\nWJdhCvtTXnm9WPfGuvhrBKHLjI6jaRZhSnEvDFxI9Ppiwnv/EJFCQDul1AxAUtNXy3gCzwbSZ3kf\nVnRZQXmP8kbHSROf/D6s6rqKt1e8zfbz242Oo2lmZ64Lqt8C6Z5P9/f3/+drPz8/m9+jMDMKvRlK\npz868XvH36lRqIbRcdKleqHqzG8/n46LO7LjrR2UzlPa6EialqzAwEACAwNNapvirZAiUgvwV0o1\nT3g9FFBKqUmJ2px++iWQD4gC+hI/RZNs30TH0LdC2rjbD2/z4pwXGVZ3GG9VfcvoOGYzY98MpuyZ\nwq7eu8idLbfRcTTNZOm6z11EHIEw4i+KXgH2Al2UUiHPaT8XWKGUWpqavrq427bo2Gia/9acKp5V\n+KrZV0bHMbtBawdx5PoR1nZbqzcA0TKMdC0/oJSKBQYA64FjwCKlVIiI9BORvkl1SalvGr4HzWAf\nrP2ArE5Z+aLJF0ZHsYgvm36Jq7Mrg9YNMjqKppmFfkJVS9HcQ3P5YucX7OmzhxwuOYyOYzF3H92l\nxuwajPEbQ9dKXY2Oo2kp0ssPaGl25NoRGv/amC09t1DBo4LRcSwu6GoQL897mcA3A6mYv6LRcTQt\nWXpVSC1NIh9H8tofr/FNs28yRWEH8C3gy5dNvqTD4g5EPo40Oo6mpZk+c9eSpJSiy5Iu5HDJwazW\ns4yOY3X9VvTjzuM7LOqwCJEkT4w0zXD6zF1LtZkHZhJ6M5Qpzaf857N79yA62oBQVjSlxRRCboTw\nS9AvRkfRtDTRxV37j7CbYYzaPIpFHReRzTnbvz7btAlKlYJixWDkSDh3zqCQFpbVKSu/tf+NjwI+\n4uTtk0bH0bRU08Vd+5fo2Gje+OsNxviNoVy+cvFvKgU3bzLnm0i6dI5j8ey7bFxyh8hIqFYNPvjA\n2MyWUsmzEqPqj+KNpW8QHWvnv6podkfPuWv/MmrTKA5cOcCqrqv+mWtW0TF83GAPy8LKsbLrArzy\n3gJHRyhShMiHTtQY24bx7Q/yWssoeOUVg78D81JK0eK3FrxY+EXGNBxjdBxN+xc9566ZZOeFncw+\nOJsf2/z4r4uI69YLK094sbvPnPjCDhAbC+fOkf36Kea1XsyAX1/gSuhdg5Jbjogwt+1cZh6Yye6L\nu42Oo2km08VdA+BB9AN6/NWD6a2mUzB7wX999vW3wtA628mT7WGSfWsWvkT/6vvp/Ut97PGXr4LZ\nC/J9i+/ptawXj2IeGR1H00yii7sGwMhNI3mxyIv/2dj6yBE4dlzoUik4+f71t3IjMiszZ1oypXFe\nq/gaPvl98A/0NzqKpplEz7lr7L64m1d/f5Wj7xwln2u+f33Wsyd4l41lWPS4FI8T6liRulNe49gx\n8PS0UFgDXY+6TuUZlVneZTk1C9c0Oo6m6Tl37fkexzzmrWVvMaX5lP8U9suXYfly6Pe2af/olit4\nl06dYNo0SyQ1Xn63/Hzb/Ft6LevF45jHRsfRtGTp4p7Jjds6jnL5yvFahdf+9+batbBiBVM/OEG3\nl86QZ9cqk4/3f/8HP/wAD5Oens/wXq/4Ot55vRm3NeXfZDTNSLq4Z2JBV4OYdWAW01pO+/cj9kFB\nRO06wuzVhfi/Usvh0CGTj+nlBbVqwbx5FghsA0SEqS2nMvPATI5dP2Z0HE17Ll3cM6k4FUe/lf34\nrPFn/7k7BmD+kcrUK3ae0nkiUn3swYPh668hLs4cSW1PoeyFGOM3hv6r+hOn7PSb1DI8XdwzqdkH\nZuPk4PTc7fKWh3vTxedomo7doAG4ucGaNelJaNv6Ve9HdGw0cw/NNTqKpiXJpOIuIs1FJFREwkXk\nPxthi0gbEQkSkUMisldE6iT67Gziz8wZXkuba/evMWrzKGa0moGD/PevwOMYR7adK07jUmdSd+Cr\nV2HqVGTaVAb7rOfrwRdg9WozpbYtjg6OzHxlJsM3Ded61HWj42jaf6RY3EXEAZgKNAMqAl1EpNwz\nzTYopXyVUlWB3sCcRJ/FAX5KqapKKX3/mA34MOBDelXpRSXPSkl+vuNsYSp43HjuQ0vPFRMDN2/C\nzZu8VmwPYVdycuiIoxkS2ybfAr686fsmg9cNNjqKpv2HKWfuNYETSqlzSqloYBHQNnEDpdSDRC/d\niS/oT4mJ42hWsOnMJrad28anDT59bpv1J0rStPSpdI2TxTGW/jX2M3tT6XQdx9aNbjCabee3seXs\nFqOjaNq/mFJ0CwMXEr2+mPDev4hIOxEJAVYAiSdyFRAgIvtE5O30hNXS50nsE95b/R5Tmk/BLYvb\nc9utP1GCZqXTv8ztG5WPsHh3UZ48SfehbJZbFje+bvo1A9YMICYuxug4mvYPJ3MdSCn1N/C3iNQF\nxgNNEj6qo5S6IiIexBf5EKXU9qSO4e/v/8/Xfn5++Pn5mSueBny/53tK5CpBG+82z21z7RqcuZ2L\nmoUvpXu8ErnuUKHwXdauzU+b5w+Z4bUv354fDvzA9H3Tef/F942Oo9mxwMBAAgMDTWqb4vIDIlIL\n8FdKNU94PRRQSqlJyfQ5BbyglLr9zPujgUil1NdJ9NHLD1jQlcgrVJpRiZ29d+KV1+u57X77DZZM\nDGNpx4VmGXfWpVZsuPsCixeb5XA2K+RGCPV/rk/wO8F4utvh2guaTUrv8gP7gDIiUlxEsgCdgeXP\nDFA60dfVgCxKqdsi4ioi7gnvuwFNgeRXoNIsYujGofSp1ifZwg6wbh00LXvWbOO+VusC69bBXftb\nDfhfynuUp6dvT4ZuHGp0FE0DTCjuSqlYYACwHjgGLFJKhYhIPxHpm9Csg4gEi8hB4HugU8L7nsB2\nETkE7AZWKKXWm/270JK188JONp7eyMj6I5NtpxSsXw9Ny6byFshk5HaPpnFjWLLEbIe0WaMajGL9\nqfXsurDL6CiapleFtHexcbHUnFOTIbWH0LVS12TbBgVBx45wos8k8y0OU748S51fZ+rU+P1X7d28\noHlM3TeVXb13JfkMgaaZk14VMhP7NehXsjplpYtPlxTbrl8PTZuaP0OrVvH/cFy4kHLbjK5b5W4o\npVhwdIHRUbRMThd3Oxb5OJIRm0bwbbNv/70w2HOsWwfNmpk/h4tL/G8EC81zjdamOYgD3zb/lmEb\nhxH1JMroOFompou7HZu0YxKNSzXmhcIvpNg2Nhb27IF69SyTpVu3+DtxMoOXir5E3WJ1mbxzstFR\ntEzMbPe5a7bl/N3zzNg/g6D+Qck3/PlnuHWL0Kt5KZi1Dbl/mmuRxdjr1oXr1+HECShb1uyHtzmT\nXp5E1ZlV6V21N0VzFjU6jpYJ6TN3OzV0w1AGvDCAIjmKJN8wKgoiI9l/Iic1ClyCyEjzBrl6Fdas\nwWHdGl6tfo6lE8PgmP2vg14sZzHee+E9hm0cZnQULZPSxd0O7b64m63ntvJxnY9N7rP/ciFqFLps\n/jAREfHzPXv20D7PFpasc4ezZ80/jg36uM7HbD67mf2X9xsdRcuEdHG3M0ophqwfwriG45JdP+ZZ\n+69YqLgn0qD4WU5H5Ob89awWHcdWuGdxx7+BPx+u/xB9m69mbbq425m/Q/8m8nEkPXx7mNwnOtaB\nI9c8qVrgigWTgbNjHG28w/hrVwGLjmNLelXtxc0HN1kZvtLoKFomo4u7HYmOjeaTDZ8wuclkHB1M\nX0f92I38FM95l+wull++sX35EJbszDzF3cnBiS+afMHHGz7Wq0ZqVqWLux2ZfXA2xXMVp2np1D2J\nZLH59iS8XOo0R87m4No1qwxnE1qUaUHh7IWZc3BOyo01zUx0cbcT9x7fY+yWsUxuMtmkB5YSs2Zx\nz+oUQ4vq1/n7b6sMZxNEhMlNJjNmyxgiH5v5biRNew5d3O3E5B2TaVamGVUKVEl1X2sWd4AOL11l\n6VKrDWcTqhasysulXuarXV8ZHUXLJHRxtwNXIq8wff90xjUcl+q+j6MdOH7DgyoFrlogWdKaV7/B\n7t1w+3bKbe3JuIbj+H7v91y7n4nmpDTD6OJuB8ZuGUuvKr0olrNYqvsevZyXMnlu4+ocbYFkSXPP\nFkvDhrBqldWGtAklcpWgR+UejNua+n+ENS21dHHP4MJvhfNnyJ8Mq5u2JyH3n8/PC4XSv6VearVt\nC8uWWX1Yw42oP4JFwYs4eTv9e9RqWnJ0cc/gRm4ayeBag8nrmjdN/fefz2/V+fanXnkFAgLg0SOr\nD22ofK75+L9a/8eozaOMjqLZOZOKu4g0F5FQEQkXkU+S+LyNiASJyCER2SsidUztq6Xdvkv72HFh\nBx/U+iDNxzCquHt4gK8vbNxo9aENN6jWILac3cKByweMjqLZsRSLu4g4AFOBZkBFoIuIlHum2Qal\nlK9SqirQG5iTir5aGiil+GTDJ4xuMBpXZ9c0HePhQwi/npPKnsZc4MusUzNuWdwYVX+UXlRMsyhT\nztxrAieUUueUUtHAIqBt4gZKqQeJXroDcab21dIm4HQAlyIv8VbVt9J8jKAgKOd5BxenWDMmM13b\ntrB8OcTFpdzW3vSu1ptTEafYfGaz0VE0O2VKcS8MJN4g7WLCe/8iIu1EJARYAbyVmr5a6iilGL5x\nOOMbjsfJIe1L8gcHQ+XCN82YzESxsfD4MWWKPiZf3jj2bHsCMZnr0fwsjlkY13AcwzYO04uKaRZh\nts06lFJ/A3+LSF1gPNAktcfw9/f/52s/Pz/8/PzMFc+uLAlZQpyKo0OFDuk6TnAwVCwYYaZUqXDw\nYPwfoG2eRiwb4UDtz+Kgfn3rZzFQZ5/OTNoxiWVhy2hXrp3RcbQMIDAwkMDAQJPaSkpnDSJSC/BX\nSjVPeD0UUEqpScn0OQW8AHiZ2ldElD6DSVlMXAyVZlTim2bf0LxM83Qdq0kTGFR2JS3zG7fe+L5L\nhej+V3tClx7PdMUdYGX4Sj7Z8AlH+h9J1WJvmgbxS1sopZJcb8SUaZl9QBkRKS4iWYDOwPJnBiid\n6OtqQBal1G1T+mqpMy9oHvnd8tOsdPp3sj52DHwK3jJDqrSrXugKkU+yEHoum6E5jNKqbCtyZ83N\nb0czyQazmtWkWNyVUrHAAGA9cAxYpJQKEZF+ItI3oVkHEQkWkYPA90Cn5Ppa4PvIFB7HPMZ/iz+f\nNfos1YuDPev2bbh/H4rmvm+mdGnjIIo2XmEs25G2+/QzOhFhYuOJjA4czZNYyy+5rGUeKU7LWIue\nlknZlN1TCDgdwMqu6d/4Yds2+Ogj2N19Gty4YYZ0abfuZGnGHG7DzuCchuYwUovfWtDaqzXvvvCu\n0VG0DCS90zKaDYh6EsXnOz5nQqMJZjnesWNQsaJZDpVuDUueJeScK1csuxGUTRvfcDwTtk3gQfSD\nlBtrmgl0cc8gvtvzHQ2KN8C3gK9ZjnfsGPj4mOVQ6ZbFMZbmNSNYscLoJMapXqg6tYvUZtreaUZH\n0eyELu4ZwJ1Hd/h699eM8RtjtmPa0pk7QLt6tzLVBh5JGdtwLJN3Tubuo7tGR9HsgC7uGcCXO7+k\njVcbvPN5m+2YwcG2VdxbvHib7dvh3j2jkxingkcFWpRtwTe7vzE6imYHdHG3cdejrjNj/ww+bfCp\n2Y554wY8eQKFCpntkOmWwy2WOnVg7VqjkxjLv4E/3+/9npsPDHhyWLMrurjbuM+3f05Xn64Uz1Xc\nbMd8Ot+ezrspza5dOzL91EzJ3CV5veLrTNr+3GcENc0kurjbsIv3LvJL0C8MrzfcrMe1tfn2p9q0\ngTVr4n+ryMxG1BvBj4d+5EpkJr59SEs3Xdxt2Pit4+lTtQ8Fsxc0zwHDw+HAAYI336BijvNw4ED8\nur+24PZtCj46Q7mSj9iy6AqcOQORkUanMkThHIXpVaUXE7aZ57ZXLXPSDzHZqNMRp3lh9guEDwhP\n8y5L//Hrr3D6NPXn9mJ0g0AalzpjnuOa0aTtdTh/NyfTWq2OP5WvVs3oSIa4HnWd8tPKc7DvQbNO\nyWn2RT/ElAGN3TKWgTUHmq+wJ1AKjt3wwCf/dbMe11zalQvl77ByxCX99zXTyO+Wn3dqvMPYLWON\njqJlULq426CQGyGsOrGKQbUGmf3Y16LcESC/W5TZj20O3vlukSvrI/Ze0sv+D6k9hGVhywi/FW50\nFC0D0sXdBvlv8WdI7SHkzGr+tVaCr+enYv7rNnenTGIdyofw5/EKRscwXO5suRlUaxCjA0cbHUXL\ngHRxtzFBV4PYem4rA2sOtMjxj133oKKHsQuFpaRjheMsCSmPvgQD77/4PpvObOLotaNGR9EyGF3c\nbcyozaP4pM4nuGVxs8jxj93IT0UP25xvf6pS/ms4OcRxMDRzrvGeWHaX7Hz80sd8Gmi+h9i0zEEX\ndxuy5+IeDl09RP8a/S02RujNfJT3sO2nH0WgY/njLNmY2+goNuHdF95l76W97L9s3I5ZWsZjUnEX\nkeYiEioi4SLySRKfdxWRoIQ/20WkcqLPzia8f0hE9pozvL0ZtXkUo+qPIqtTVouNEXYrL955bbu4\nA3SoEMKfG3PpqRkgm3M2RtQbwajNo4yOomUgKRZ3EXEApgLNgIpAFxEp90yz00B9pZQv8Ztjz0r0\nWRzgp5SqqpSqaZ7Y9mfL2S2cijhFryq9LDZGRFQWHkQ7Uyi77T8cVL3gZZ5EC8HBRiexDb2r9ibk\nRgjbz283OoqWQZhy5l4TOKGUOqeUigYWAW0TN1BK7VZKPV2ndDeQ+D42MXGcTEspxcjNIxndYDTO\njs4WGyfsSg68896y6TtlnhKBDo3u8OefRiexDS5OLnza4FNGbhqJfthPM4UpRbcwcCHR64v8u3g/\nqw+wJtFrBQSIyD4ReTv1Ee3f+lPrufngJt0qdbPoOGFXc+Cdz/anZJ7q+PIdliwxOoXt6OHbg8uR\nl9l4ZqPRUbQMwKxn1CLSEOgFJJ6Xr6OUqga0BN4TkbrmHDOje3rWPsZvDI4OjhYdK+xKTrzz3rLo\nGOb0ok8UEREQGmp0Etvg5ODEGL8xjNg0Qp+9aylyMqHNJaBYotdFEt77l4SLqLOA5kqpiKfvK6Wu\nJPzvDRH5i/hpniQnDv39/f/52s/PDz8/PxPiZWzLwpYRHRtNxwodLT5W2NUcvF4i41RKBwfo2BF+\n/x1G6+d4AHjd53Umbp/IyvCVtPZubXQczcoCAwMJDAw0qW2KC4eJiCMQBjQGrgB7gS5KqZBEbYoB\nG4HuSqndid53BRyUUvdFxA1YD4xRSq1PYpxMt3BYbFwsVWZWYWLjibzi9YrFx6tY+A4LWi/Et8A1\ni49lFm3asCe6Gm++CSEhtrf+vFGWhS7j08BPOdTvEA6iL2dlZulaOEwpFQsMIL4wHwMWKaVCRKSf\niPRNaDYKyANMf+aWR09gu4gcIv5C64qkCntm9fux33HP4k6rsq0sPlZsLJy+4U7ZvLctPpbZnDxJ\nzegdREc+5PDPh2HHDoiISLmfnWvj3QYXRxf+PK6vNmvPp5f8NUh0bDQVpldg5iszaVSykcXHO3UK\nGte6z9n3vrT4WOY2clMjHsc4MrlpAHTrBmXLGh3JcOtPref9Ne8T/G4wTg6mzK5q9kgv+WuDfgn6\nhWI5i1mlsEP8RUnvAhlz9+kuPkdZdMwn0y8DnFiTUk3wdPdk/pH5RkfRbJQu7gZ4HPOYcVvHMb7h\neKuNGRYG3gXvptzQBlXMf4M82R6y/XyxlBtnEiLC+IbjGbNlDI9jHhsdR7NBurgbYOaBmVTKX4na\nRWtbbcywMChXMGOeuQN08QlmwdFKRsewKfWK16N8vvLMOTjH6CiaDdLF3cruP7nPZ9s+Y0Ij6+6P\nGT8tkzHP3AE6+wSzJKR8pt88+1njG41nwrYJRD2xzc1XNOPo4m5l3+35joYlG+JbwNeq48ZPy2Tc\nM/cSue7gnfcWATtcjY5iU6oVrEadYnWYuneq0VE0G6OLuxVFPIzgm93fMMZvjFXHvXMHoqKgcO4H\nVh3X3Lr4HGXByuxGx7A5Y/3G8tWur7j7KOP+ZqaZny7uVjR552TaebfDK6+XVccNCwMvr4z/ENDr\nPsdYFejGXV3D/qW8R3laebXiy50Z7zZXzXJ0cbeSq/evMvPATD5tYP0ddcLCwNvb6sOaXT7XBzSp\n84BFi4xOYntGNxjN9P3TuR5l27tsadaji7uVTNg6ge6Vu1M0Z1Grjx0WBuWeXYE/g3qrwz1+/NHo\nFLanRK4SdKvUjc+2fWZ0FM1G6OJuBWcizrAgeAEj6o0wZHx7OXMHaFr3AZcvw1G9X/R/jKg3gnlH\n5nHuzjmjo2g2QBd3KxgdOJqBNQfi4eZhyPihofZT3B0doWdPmDvX6CS2x9Pdk3drvIv/Fn+jo2g2\nQBd3Czt67SjrTq1jcO3BhowfGxu/rozdLMdy/jw965xg/i8xPDl8HI4fh4cPjU5lMz586UNWha/i\n+I3jRkfRDKaLu4WN2DSCYXWHkcMlhyHjnzsHHh7g5mbI8Oa3bRtl9vxGhewXWTnxKCxeHH+vpwZA\nzqw5+eilj/Rm2pou7pa04/wOgq4F0b9Gf8MyhIfbz5RMYm9VPcSPh6oaHcMmDag5gD0X97Dn4h6j\no2gG0sXdQpRSDNs4jNENRpPVKathOezpYmpiHcofZ9eFoly4a8xvRLYsm3M2/P38GbpxqN6OLxPT\nxd1CVp1Yxa2Ht+jh28PQHOHh8Q8w2Ru3LNG8UfkIP+yvYXQUm9SzSk+u3r/KulPrjI6iGcSk4i4i\nzUUkVETCReSTJD7vKiJBCX+2J+ynalJfexQbF8vQDUOZ2HiicRspHDwIy5cTtv0G3rd3wfLlcOOG\nMVksZEDNvcw+WF1fT02Ck4MTnzX6jE82fEKcijM6jmaAFIu7iDgAU4FmQEWgi4g8+0jMaaC+UsoX\nGE/8Rtmm9rU7847MI3e23LT2MnAD47Nn4eBBws5mwStiT3yxj4w0Lo8FeOW9RY1Cl1m0LJvRUWxS\nu3LtcHV2ZcHRBUZH0Qxgypl7TeCEUuqcUioaWAS0TdxAKbVbKfV0xY/dQGFT+9qbh9EP+XTzp0x6\neRJi8GIuUU+cufnAlWI57Xcxlvdf3MN3P7qhp5b/S0SY9PIkRm0epTf0yIRMKe6FgQuJXl/kf8U7\nKX2ANWm358rKAAAgAElEQVTsm+FN2zeN6oWq81LRl4yOwonbeSmT5zaODvZb+ZqWPkXUQ2HHDqOT\n2Kb6xetT0aMiM/bPMDqKZmVmvaAqIg2BXkCmmFt/VsTDCCbtmMRnjWxjfY+wm3nxynvL6BgW5SCK\ngb2i+O47o5PYrs9f/pyJ2ydy55F+HiAzMeVq3yUg8eaVRRLe+5eEi6izgOZKqYjU9H3K39//n6/9\n/Pzw8/MzIZ7t+GzbZ7xa7lXKe5Q3OgoA4bfy4m3nxR3gzU4PGf11Ti5cgKLWX5fN5vnk96G1V2sm\nbpvIpCaTjI6jpUNgYCCBgYEmtZWU7oMVEUcgDGgMXAH2Al2UUiGJ2hQDNgLdlVK7U9M3UVuVke/J\nPXvnLNVnVSf4nWAKZi9odBxYupQ3/MvwcqnT9Kxy2Og0ltWvHx98XhAXF/jiC6PD2KZL9y5R+YfK\nHOp3iGI59Ubj9kJEUEoleXEvxWkZpVQsMABYDxwDFimlQkSkn4j0TWg2CsgDTBeRQyKyN7m+6f6O\nbNDwjcN5v+b7tlHYE8Sfud80OoZVDBkCc+bALfv/RSVNCucozLs13mXkppFGR9GsJMUzd2vJyGfu\n+y/vp+2itoQNCMM9i7vRcQBQS5aSq1tLznwwhTzZ7PxGcB8fcHWlz5RKFMrziLHdT0DjxuDiYnQy\nmxL5OBKvqV6s7rqaqgX10g32IF1n7lrylFJ8uP5DxviNsZnCDnDtjgvODnH2X9gBgoNh716GeS9l\n+rLC3N0aBNHRRqeyOdldsvNp/U/5KOAjvSxBJqCLezqtCF/BzQc36VWll9FR/iX8sjve+TLHlMxT\npfNE0LLsCaburWl0FJvVp1ofLkVeYvWJ1UZH0SxMF/d0iI6N5uOAj5ncZDKODo5Gx/mXsMvZ7f42\nyKQMr7eNKXtqcf++0Ulsk7OjM182+ZIPAz4kOlb/dmPPdHFPhxn7Z1AiVwlalG1hdJT/CL/snilu\ng3xWuXw3aVTyDD/86Gx0FJvVsmxLiuQowqwDs4yOolmQLu5pdPvhbcZvHc9XTb8yOkqSwi5nzzR3\nyjxrRL2tfPW9M/fuGZ3ENokIXzX9irFbxxLxMCLlDlqGpIt7Go0JHEPHCh2pmL+i0VGSFHbZPVNO\nywBU8rxOs8axTNLP6zxXZc/KtPFqw4RtE4yOolmIvhUyDcJuhlF3bl2Ov3vcsE2vkxMdDdndYrn7\n8We4OMUaHccQl7p+ROXabhw8CMWLG53GNl29fxWf6T7s7rObMnnKGB1HSwN9K6SZfRjwIZ/U+cQm\nCzvAmTNQKPejTFvYAQoXUgwcCMOHG53EdhVwL8CQ2kP4cP2HRkfRLEAX91Rad3IdITdCGFhzoNFR\nnis8HLwL29fa7Wnx0UcQGAh79FaizzWo9iCOXj9KwKkAo6NoZqaLeyo8iX3CB2s/4Jtm3+DiZLtP\nP4aGgnehTH4v4IIFuC36kfFNtzK42zXUnB+NTmSTsjpl5aumX/F/6/5P3xppZ3RxT4Wpe6dSMndJ\nXvF6xegoyQoJgfJFMvmZ++XLcOECPYpu5mFUHPNX5TY6kc1q692WQtkL8cP+H4yOopmRLu4munb/\nGhO3T+SbZt8YvsNSSkJCoHxhfR8ggKODYnbr5XwY0JSrV41OY5tEhG+bfcu4reO4+SBz3j5rj3Rx\nN9HwjcN50/dNyuWz7S1gldJn7s+qXugKvase4t130dvxPUfF/BXp7NNZrxppR3RxN8G+S/tYfXI1\no+qPMjpKiq5fBwcHyJfjidFRbMqnDbYQGgp//GF0Ets1xm8My8KWceDyAaOjaGagi3sKYuNieXf1\nu0xsPJGcWXMaHSdFISFQvjzY+MyR1WV1imHuiJO8/240N7Ych6NH4688a//InS03nzX6jPdWv0ec\nijM6jpZOurinYM7BObg4utDDt4fRUUzytLhr//Xiifl099pLv/6C+nMJrFpldCSb82aVN3EQB+Ye\nmmt0FC2dTCruItJcREJFJFxE/rP5tYh4i8hOEXkkIoOf+eysiAQl3qEpo7gRdYNRm0cxvdV0HCRj\n/Duoi3vyxjfaxKXI7EzeWcfoKDbJQRyY1nIawzcN5/bD20bH0dIhxYolIg7AVKAZUBHoIiLPXlW8\nBQwEJidxiDjATylVVSmVoRbaHrphKN0qdaOyZ2Wjo5hMF/fkuTjF8udri/lmdy02ntB7iSalasGq\nvFbhNUZsHGF0FC0dTDkdrQmcUEqdU0pFA4uAtokbKKVuKqUOADFJ9BcTx7Epuy7sYu2ptYxpOMbo\nKKkSGqqLe0qK5rzHb+2X8saCFpw/b3Qa2zSu4Tj+DvubvZcy1C/bWiKmFN3CwIVEry8mvGcqBQSI\nyD4ReTs14YwSHRvNO6veYXKTyeRwyWF0HJNFRsLt21BMn5CmqFHJMwxpcIAOHSAqyug0tid3ttxM\nbjKZfiv7EROX1DmbZuucrDBGHaXUFRHxIL7IhyiltifV0N/f/5+v/fz88PPzs0K8/5qyZwr53fLT\nxaeLIeOnVWgoeHnF3wqppWxIg/0cP96Adu1gxQrImtXoRLalW6Vu/Hz4Z77b8x2Daw9OuYNmcYGB\ngQQGBprUNsUlf0WkFuCvlGqe8HoooJRS/1ktW0RGA5FKqa+fc6znfm4rS/6evXOWGrNqsKfPHkrn\nKW10nFT59VdYuxYWLACWLoUjR4yOZNuyZyf2/4bQtSs8fAhLloCz3sDpX8JvhfPSjy9xsN9BiuXU\nvxLamvQu+bsPKCMixUUkC9AZWJ7ceIkGdhUR94Sv3YCmQLDJya1MKcW7q95lSO0hGauwKwVKEXJc\nUb6c0o9hmurxYxy3bGJ+781w8wbdmlwjZtc+o1PZFK+8Xrz/4vsMXGO7q6BqSTNpsw4RaQ5MIf4f\ngx+VUp+LSD/iz+BniYgnsB/ITvzdMfeBCoAH8Bfx8+5OwG9Kqc+fM4bhZ+6Ljy1m7JaxHOx3kCyO\nWQzNkioBAbBjB6/+/jrdKh2lY4XjRifKcB7FONFmYRdcszswf3tJ3N2NTmQ7Hsc8xvcHXyY2nsir\n5V81Oo6WSHJn7ibNuSul1gLez7w3M9HX14CiSXS9D1QxPapxIh5GMGjdIP547Y+MVdgTCbnhQbl8\neuGntMjqFMOKLgt4Z/Nr1K0Ly5dnrAvTT57A1asQEwNxCQ+XFisGWczwV9nFyYVZrWfRZUkX/Er4\nkTubXmEzI7DGBdUMYcj6IbTzbsdLRV8yOkqaPIl15OydXJTNkzn3TTUHF6dYfnzlb7463oLaVcqy\nZPBOanndhoIFoV49o+P9IyYG9u2DjRvh0CE4fjx+9618+eKvGThKLLH3orh2LxsVCkRQregNWtSP\nos34mjg6pm3M+sXr09a7LR8FfMScNnPM+w1pFqGLO7D+1Ho2ndnE0XeOGh0lzU7ezkOxnHcz9dZ6\n5iCPH/Fh6b/wbuZFm4lt6VPtIJ/2vkRWg2v79evxv00sn3uTLUG5KFngEY2rRfC6zz0qtH5I2SIP\nccmSMK15/z7s3UvUE2eCrhXgwOWCfL64Jh/9AYMHQ8+e4Oqa+gyfv/w5PtN92Hh6I41LNTbr96eZ\nX6bfIDvycSSVZlRiVutZNC3d1Orjm0VAAEtm3+bXI74s67zI6DR240qkOwPXtOTonSLMXpSD+vWt\nO/7Fr35nyc6CLAkqTdClfDQrd552pY7wcqnT5HdL3c35yiM/O9yb8dWv+dhz1JVfxl+kSUtnKFIk\nVcdZfWI1A9cM5Ej/I7hlcUtVX838kptzz/TFfcDqATyIfsBPbX+y+thmExDAhIkORD5x4fOXNxid\nxu78dbcRA/+oT7VqMGwY1K5tgUGioiA2lpOnhGWrnVmy3Imw4GjaeIXSvlwITUqfJquTeR4m2nSm\nJN3/epW3ml1m9MJyOKXy9/fuf3XHw9WDr5slecezZkW6uD/HtnPb6LykM8HvBGfsi0QBAbwxxJMm\npU7xZpUgo9PYHy8vHr7alZ9/hsmToWhRGDAAmjeH7NnTcLzdu+HKFQAePnZgV3geAjY7syy4NBGP\nstLaK5z25UNoVPIMWRwtM8127b4bbwS8yZNc+fnzT/DwML3vzQc3qTyjMotfW0zdYnUtkk8zjS7u\nSbj/5D6+P/jyTbNvaOPdxmrjWkRAANV6+PDDKyupWfiS0Wnsj5MTT++NjIkVFh/24pfg6uw6mY+6\nLzymVaNHVKkYTYVqWcld/PnLVURHQ1gYHJ25kyOHYtl5sSgHLheisuc1GpU8QxvvMGoUuoyDWOe/\ng1jvCgzf3Yb1m53YvPw+uXIqcHExaUJ+edhyBq8bTFD/ID09YyBd3JPQf2V/Hsc+Zm7bjL9udfSa\nDeRs04CbH3+Bq7Pewd5a7j5yYd2pMqw7WZpjN/JzPKIAbjmc8PCIX8ogWzZQcXHcvAk3bwl37kDJ\nEopKuS5S2e0UNQtfom6x87hnMW7XLKXgg7UtOHilIOvemIdboxfh5ZdN6vvm32+SPUt2pracauGU\n2vPo4v6MNSfW8M6qdwjqH5QhdldKybHZO2k/zIuwAfo/MiOpsl5czOVDRKQTDx8LDx85wOnT5Is8\ng4drFHldH+LkYHs7HMUpodeytly9787y787h0tK0O2HuPLpD5RmV+antT7xcyrR/EDTz0sU9kdsP\nb1N5RmXmvTqPhiUbWnw8a1gw7Ch/LXfgj9f0BqFa2sTEOdDpj9dwK5yLXwMKmrxN47qT6+i7si9H\n+h+xixOljCa9a8vYDaUU761+jw7lO9hNYQc4csYdX89rRsfQMjAnhzjmt1/KwZPZmTfP9H7NyjSj\nVdlWvLf6PcuF09IkUxX3+UfmE3Q1iIkvTzQ6ilkFncpOZV3ctXRydY5m4fBghgyBkydN7/dl0y85\neOUgvx35zXLhtFTLNMX91O1TDF4/mIUdFuLqnIbH82xY/Jn7VaNjaHagcqn7fPopdOkSv16NKVyd\nXVnQYQGD1g3iTMQZywbUTJYpint0bDRdl3ZlZL2R+BbwNTqOWd28CVGPHCmW867RUTQ7MWAAeHrC\np5+a3qdKgSp8UucTuv/VXe/cZCMyRXH3D/Qnb7a8vP/i+0ZHMbsjR6ByyfsmXwDTtGTt2YN8OZm5\n1afyy/T77B04Dw4cMKnroNqDyOacjQlbJ1g4pGYKuy/um85sYu7hucxtOxexwwp45AhULhVpdAzN\nXkRHQ1QUHnKTzxsFMGBpI+Iem/bshIM48Gu7X5l5YCabzmyycFAtJXZd3K9EXuGNpW/w66u/4unu\naXQciwgKij9z1zRz6+57BGfHWH5aafraBAWzF+SXdr/Q/a/uXL2vrwMZyaTiLiLNRSRURMJF5JMk\nPvcWkZ0i8khEBqemr6XExMXQZUkX+lbva9cPWBw5Ar76zF2zAAdRTGu5mhE/FOX2bdP7NSndhLeq\nvMUbS98gNk4vQW2UFIu7iDgAU4FmQEWgi4iUe6bZLWAgMDkNfS3CP9AfJwcnRtUfZY3hDBETAyEh\n4FNCn7lrllGlwFU6NrzNyJGp6zfabzQxcTFM2Kbn341iypl7TeCEUuqcUioaWAS0TdxAKXVTKXUA\nePYyeYp9LWHtybX8fPhnfmv/G44Oadx6JgMID49fjtstm+090q7Zj/H9LrB0afyuT6ZycnBiQYcF\n/LD/B9adXGe5cNpzmVLcCwMXEr2+mPCeKdLTN01OR5zmzb/fZEGHBXY7z/7UkSNQubLRKTR7lztH\nLP7+8EkqJ1ULZS/Ewg4L6fF3D05HnLZINu35bGqbPX9//3++9vPzw8/PL1X9o55E0W5RO0bWG0n9\n4lbeNscAQUG6uGvW0bs3fPklbN4MDVOxckeDEg0YXnc47X9vz87eO+3uAUJrCwwMJDAw0KS2KS4c\nJiK1AH+lVPOE10MBpZSalETb0UCkUurrNPRN18JhSik6L+lMNqdsdnvb47NatYK+faGtawDs2GF0\nHM1eeXtD8eL8ti4fU5cUZOfMo4i3V/yO3CZQStH9r+4oFPNfnZ8p/tu0lvQuHLYPKCMixUUkC9AZ\nWJ7ceOnom2Zf7PiC0xGn+eGVHzLNXx595q5ZRVgYrF9PFxYSdeMBK6eehWumr2UkIsxqPYvjN47z\n9S69NZ+1pFjclVKxwABgPXAMWKSUChGRfiLSF0BEPEXkAjAIGCEi50XE/Xl9zf1NrAhbwZQ9U1ja\naSlZnbKa+/A26dYtiIyEEiWMTqJlFg6iGN9oEyM2NSYuldfwXZ1dWdZ5GV/t+ooVYSssE1D7lwy/\nnvvhq4dpMq8JK7us5MUiL1ogmY2JiIAbN1i71ZXPZ+chcN7F+NsYQsz+b6am/YdS8NJPvRk4UOj6\ncZFU9999cTetF7ZmY4+NVPbUv3aml92u53458jJtFrZhWstpmaOwA4SGwoIF7J5/ktrOB2DBAl3Y\nNasRgQmNNjF6mgcxaVgfrFaRWnzX/DvaLGzDtft6mWpLyrDFPepJFG0WtqF/jf50qtjJ6DhWt+ti\nEWoXvWh0DC0TaljiDAXyxbBoUdr6d6nUhTd936TNojZEPYkybzjtHxmyuEfHRtN5SWcqeVZiWN1h\nRsexujgl7LlYhFpFdHHXrE8EPu1/gwkTIDaNqwv4+/lTPl95Oi/prJcItpAMV9yVUvRd2ZfYuFhm\nvTIr09wZk1jIjXzkc31Afjd91qMZ4+XaUeTMCUuWpK2/iDC79WyiY6Ppv7I/tnLtz55kuOI+fONw\njt84zh+v/YGzo7PRcQyx62JRahe9kHJDTbMQERg1CsaPJ9V3zjzl7OjMn53+5PDVw/gH+ps1n5bB\nivuU3VP4K/QvVnVdhVsWN6PjGGbXhSLU1lMympGOH6el+1acH0eyfOIx2LoV7qZ+NzD3LO6s6rqK\nBcEL+G7PdxYImnllmOI+5+Acvt79NeveWEc+V9OejLNXuy4WpXYRfeauGejYMWTzJkZVXcW4ablR\nGzfBnTtpOpSnuycbum/gq11f8ePBH80cNPPKEMV9XtA8/AP92dB9A8VzFTc6jqEi7jly4V4OKnle\nNzqKptHGO4zoWEfWnCybruMUz1WcgO4BfBr4KQuPLjRTuszN5ov74mOL+XjDx6zvvp6yedP3F8ge\n7DmenRqFLuPkoJf51YznIIoR9bYyfmt90ntN1CuvF+veWMegdYP4K+Qv8wTMxGy6uP95/E/eX/M+\na7utpYJHBaPj2ITdwe7UKqzn2zXb0bHCcW49zEbgLpd0H8snvw+ru63mnVXv8OfxP82QLvOy2eI+\n/8h8Bq4ZyLo31uFbwNfoODZjV3B2/fCSZlMcHRTD6m5nwrScZjletYLVWPfGOgauGainaNLBJov7\nnINzGLphKBt7bNSFPZG4ONhz3F0/vKTZnG6VjnDyrDO7d5vneL4FfAnoHsCQ9UP4NehX8xw0k7G5\n4v7dnu8Yt3Ucm9/crKdinhESAvlyxuiHlzSb4+wYx8f97zLBjFum+uT3YWOPjYzYNIIpu6eY78CZ\nhE0V96EbhjJ933S29tyqL54mYdcuqO0TaXQMTUvSW6/d58ABOHzYfMcs71Ge7b22M2P/DIZvHK6f\nZE0FmyruW85tYcdbOzL97Y7Ps3kz1PO9Z3QMTUtS1qzw4YfxT62aU/Fcxdn+1nY2nN5An+V99Fo0\nJrKp9dyjnkTpPRafIyYGPD0h6KcDFDmkNzvQbFDjxkS55adMs9KsnX0B33KPoVQpyJLFLIe//+Q+\nHRd3BOD3jr+TM6t5LuBmZOlez11EmotIqIiEi0iSe6CLyHcickJEDotI1UTvnxWRIBE5JCJ7kxtH\nF/bn270bihaFIvmfGB1F05K2cSNuyxfyUbWNjBn2CBYtgvv3zXZ49yzurOy6ktK5S/PSTy9xJuKM\n2Y5tj1Is7iLiAEwFmgEVgS4iUu6ZNi2A0kqpskA/YEaij+MAP6VUVaVUTbMlz2RWrYJXXjE6haal\nrH+N/ey6WITDVwuY/dhODk5MazWN/tX789JPL7HjvN4Y/nlMOXOvCZxQSp1TSkUDi4C2z7RpC/wK\noJTaA+QUEc+Ez8TEcbRkrFwJrVoZnULTUubqHM3HL+1gzJYGFhtj4IsD+anNT7z6+6tM3TtVX2hN\ngilFtzCQeJWqiwnvJdfmUqI2CggQkX0i8nZag2Zm58/D1atQU//eo2UQ/WvsZ8/FIhw64mixMVqU\nbcGu3ruYc3AOPf7uwYPoBxYbKyNyssIYdZRSV0TEg/giH6KU2p5UQ39//3++9vPzw8/PzwrxbN+q\nVdC8OTha7r8TTTOrbM4xfFJnO/6TGrPMz3LjlM5Tmp29d9J/ZX9qzanF7x1/p7xHecsNaLDAwEAC\nAwNNapvi3TIiUgvwV0o1T3g9FFBKqUmJ2vwAbFZK/Z7wOhRooJS69syxRgORSqmvkxhH6V+tnnH+\nPOzZQ6uJdenR4Cyvv3QRbt2KP43XNBv3MNoJr5+H8fsfjrz0kmXHUkox5+Achm8azoRGE3i72tuZ\nYpe29N4tsw8oIyLFRSQL0BlY/kyb5UCPhMFqAXeUUtdExFVE3BPedwOaAsFp/D4yn3v3eHA4nG3H\n8tDMZQscO6YLu5ZhZHOOYcKIhwwZQrpXjEyJiPB29bfZ1msbM/bPoOMfHbn14JZlB7VxKRZ3pVQs\nMABYDxwDFimlQkSkn4j0TWizGjgjIieBmcC7Cd09ge0icgjYDaxQSq23wPdhtzafKUHVglfIlfWR\n0VE0LdXe6PSER4/gjz+sM165fOXY3Xs3JXOVpNKMSpl66WCbeojJVrLYjOBg3un5kFK5I/iozk6j\n02ha6g0YwKbDeejTVwg5pnBxIX4DVitMmew4v4Ney3pRrWA1vm/xPR5uHhYf09rS/RCTZozYWFgR\n7k0rrxNGR9G0tJk6lUbbx1LROZxp7QJg7FiIjrbK0HWK1SGofxBFcxTFZ4YPPx78kTiVeTa50cXd\nhq3amp0iOe5RweOG0VE0LV2+eDmAidvrcutBNquOm805G5ObTmZNtzXMPjibuj/V5fBVM65sZsN0\ncbdhUxfm5b0Xkl2xQdMyhPIeN+lcMZiPA5oYMn61gtXY2Xsnb1V9i2bzm9FvRT+u3rfvmxN0cbdR\nYWEQFJaV1yoeNzqKppnFhMab2HCmFAEbjLlF0UEc6FOtD6HvhZLDJQcVp1dk7JaxRD2xz/0RdHG3\nUdOnQ58Ot8nqpJc31exDDpfHzHxlJX3fdTLnemKpljtbbiY3ncz+t/cTcjOE0t+V5utdX9vdE676\nbhkbFBkJxYtD0OIwim7Xe0hq9qXn6VHkyOXId98ZnSTekWtHGLNlDDsv7OTjlz7m7epv457F3ehY\nJtF3y2Qw8+dDw4ZQtIB17irQNGv6+otYliyB7UkuQmJ9lT0rs6TTEtZ2W8vOizspOaUkIzaOyPBz\n8rq42xilYOpUGDDA6CSaZhl5Vv7KtFc30OPVSG5O+Q3mzTM6EhC/Kfcfr/3B7t67ufPoDhWmVaDn\n3z3Ze2lvhlx1Uk/L2IoLF2DhQlaHluKjVX4ED/4JiY2BJ3pzDs0+fRLwMnsvF2Z978U4j0xyDyBD\n3Xxwk7mH5jJj/wzyZMtDv+r96FSxk03tAJXctIwu7rbi7FkiZy6g8g/v8EOrlTQrc8roRJpmUbFx\nQptFXSiZN5Kpu6obHee54lQca0+uZc7BOWw6s4lWXq3o6duThiUb4uRgjYV1n08X94zg7FkGdLrO\ng2hnfmq7zOg0mmYVdx+58OKPfRnyWV7ezgC7PdyIusHC4IX8GvQrF+5doEP5DnSq2Il6xerh6GD9\nNbl1cc8Atv5+ha5vu3H0nenkzqYXCdMyj/DIgtSb349p06BjR6PTmO7U7VP8cfwPfj/2OxfvXaRl\n2Za09mpNs9LNyO6S3SoZdHG3cQ8egG/FaL6q9SdtvMOMjqNp1uXiwuEWw2jVCkaNgv79jQ6Ueufu\nnGNl+EpWhK9gx4UdVC1QlSalmtCkdBOqF6yOs6OzRcbVxd2GxcXF/2W+f+0+C6p+aXQcTbM+R0do\n2JBTl7LSdHBF3mx+nVHv3ERq2O48fHKinkSx7fw2Ak4FEHA6gLN3zlKrSC3qFatH3WJ1qVGohtnO\n7HVxt1FPnsBbb8GZM7By+nly//WT0ZE0zVBX77vTfP4bVChylylrvPGwg1V6bz24xY4LO9h2bhs7\nL+4k6GoQJXKVoGbhmlQrWI0qBarg6+mbpoKf7uIuIs2Bb4m/L/7HxFvsJWrzHdACiAJ6KqUOm9o3\noV3mKe4REdybtYgOP7XC3eUJC7qvJZs8wtBnsjXNRkQ9cWb0nhbMC67G5MnQvbtVln+3mujYaI5e\nP8reS3s5dOUQh68dJvh6MJ5unlTMX5GKHhWp4FEBr7xeeOX1Ik+2PM89VrqKu4g4AOFAY+Ay8dvu\ndVZKhSZq0wIYoJRqJSIvAlOUUrVM6ZvoGIYX98DAQItvyq0UbFt5l/fffkjtIheY2nI1jg62949a\n4Nmz+JUoYXQMm6B/Fv9jtZ+FqysH8jbh7UllyeEWy+DXL9GyWSxOVStZfmwTmbNexMTFcDriNMeu\nH+PYjWMcv3GcE7dPEH4rHCcHJ0rlLkWp3KUomaskxXMWp2jOohTLWQzfAr7PLe6m3KRZEzihlDoH\nICKLgLZA4gLdFvgVQCm1R0RyiognUNKEvjbDksVdKVi7FiZMgKuX3RlZazNv+h622TMSXdD+R/8s\n/sdqP4sHD6j+YBl7Ozuw4GglJk6rwXtf5OHt/4O2baFSJXAw+Pl6c9YLJwenf87UXy3/6j/vK6W4\nHnWd0xGnOXPnDGciznD46mGWhy/nwt0LyR/ThHELA4mPcpH4gp9Sm8Im9rU7jx/DlSvxD53u2xXD\njkUX2HHKk/zZHzK8xWE6dj6B082MvW6FplmDk0McPXyD6OEbxOF7pZizrwGdZuXl+l0X6laI4IVq\nsZSuW5BSpaBECcibF7JkMTq1+YgInu6eeLp7Urto7f9+/t7zzw4t9XhVms5HW7c2d4zUCQuDAwfi\nv+kEkXEAAAOaSURBVH52hkjFxaFU/N0tSsVvgRfzRBEbC48ex9/O+OChcO9OHJFRDhT0iKGQRwzV\nykXRoeppvu64k2J57iecqQsULGjtby913N1tP6O16J/F/xj4s6hS8CFTvdfCK3D1bja2nSzIoSuF\nWfmzE6cvZ+XsFRdu33MiSxbInVtwd4/fs9UlC7hkBUdHwdEx/uacp9u4ioP8Z0vX5339rMT1whaZ\nMudeC/BXSjVPeD0UUIkvjIrID8BmpdTvCa9DgQbET8sk2zfRMWxv4lnTNM3GpWfOfR9QRkSKA1eA\nzkCXZ9osB94Dfk/4x+COUuqaiNw0oW+yATVN07TUS7G4K6ViRWQAsJ7/3c4YIiL94j9Ws5RSq0Wk\npYicJP5WyF7J9bXYd6NpmqYBNvQQk6ZpmmY+erOO5xCRISISJyLPf4LAzonIFyISIiKHRWSJiOQw\nOpM1iUhzEQkVkXARsb0Fx61ERIqIyCYROSYiR0XkfaMzGU1EHETkoIgsNzrL8+jingQRKQI0Ac4Z\nncVg64GKSqkqwAlgmMF5rCbhAbypQDOgItBFRMoZm8owMcBgpVRFoDbwXib+WTz1AXDc6BDJ0cU9\nad8AHxkdwmhKqQ1KqbiEl7uBIkbmsbJ/Ht5TSkUDTx/Ay3SUUlefLieilLoPhBD/DEumlHDy1xKY\nY3SW5Oji/gwRaQNcUEodNTqLjXkLWGN0CCt63oN5mZqIlACqAHuMTWKopyd/Nn3B0tg9ogwiIgGA\nZ+K3iP8/aiQwnPgpmcSf2a1kfhYjlFIrEtqMAKKVUgsMiKjZCBFxB/4EPkg4g890RKQVcE0pdVhE\n/LDh+pApi7tSqklS74uID1ACCBIRIX4a4oCI1FT/394do0QMBFAY/h/aCHoD8RC2CwpbicX2Iohn\n0HNYCR5gbbyDN7DZysLKSgRRPIE8i2RBUFPOhMn7uqR6hPCYCTMT+61gxGL+exZrks7ppqDzIoHG\n4wXY+3G929+bJEmbdMV+a3vK/4GcAQtJx8AWsCNpafuscq5fshRygKRnYN/2Z+0sNfTHNV8BB7Y/\naucpSdIG8ER3oukr8ACcTHWfhqQl8G77onaWsZB0CFzaXtTO8pd8cx9mRjztKuAa2Abu+2VfN7UD\nlWL7C1hvwHsE7iZc7DPgFJhLWvXvwlHtXDEsI/eIiAZl5B4R0aCUe0REg1LuERENSrlHRDQo5R4R\n0aCUe0REg1LuERENSrlHRDToGxjKaKyUxm7/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7328d4ca10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xsample, ysample, accept, c, normalize = sampling(-1.0, 2.0, nsamples=1000000)\n",
    "\n",
    "I = quad(pdftarget, -100, 100)\n",
    "xmin, xmax = -5, 5\n",
    "N = 100\n",
    "x = np.linspace(xmin, xmax, N)\n",
    "target = pdftarget(x, normed=I[0])\n",
    "propos = c*norm.pdf(x, -1, 2)\n",
    "\n",
    "plt.hist(xsample[accept], bins=50, normed=1, histtype='stepfilled', linewidth=0, color='red', alpha=0.5);\n",
    "plt.plot(x, target)\n",
    "plt.plot(x, propos)\n",
    "plt.xlim([xmin, xmax])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make an animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_samples(mu, sigma, 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_sampling/'+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",
    "    \n",
    "    x = np.linspace(xmin, xmax, 100) # draw target dist line\n",
    "    target_dist = pdftarget(x, normed=normalize)\n",
    "    propos_dist = c*norm.pdf(x, mu, sigma)\n",
    "    \n",
    "    ymax = 1.2*np.amax(propos_dist)\n",
    "    ax.set_ylim(0, ymax)\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": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAD7CAYAAACc26SuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdUFNffx/H3pSqiIIpYwG4s2BUrCrbYe4nGWGOPiUZj\nSUwsMcX4i/ExYosaa+waK/aIBo1BJTbsxhoLqNgAEdj7/IEiKsIiu8yW+zpnz2FmZ+Z+WPHLcOfO\nHSGlRFEURTFPNloHUBRFUd6eKuKKoihmTBVxRVEUM6aKuKIoihlTRVxRFMWMqSKuKIpixuwyszEh\nhBrPqCiK8haklCKl9Zl+Ji6l1Pw1btw4zTOYykt9FuqzUJ+F6X8WqVHdKYqiKGZMryIuhGgihDgj\nhDgnhBiVwvt+Qoj7QojQZ68vDR9VURRFeVWafeJCCBsgAGgA3AAOCSE2SCnPvLLpPillKyNkNDh/\nf3+tI5gM9Vm8oD6LF9Rn8YKpfxYirf4WIUQNYJyUsumz5dGAlFL+kGwbP+AzKWXLNI4l02pPURRF\neZkQApmBC5sFgGvJlq8/W/eqmkKIo0KILUKIMm+RU1EURUknQw0xPAIUlFJGCyGaAuuBd1LacPz4\n8Ulf+/v7m/yfKsrrEhIgMPDldc2aga2tNnkUxdIEBQURFBSk17b6dqeMl1I2ebb8WndKCvtcAqpI\nKe+9sl51p1iA+Hj45puX1335Jdhl6l0HirWTUjJ8+HB++umnpHUTJ06kQoUKnDx5ki+++CJd60xZ\nRrtTDgHFhRCFhBAOQGdg4ysNeCT7uhqJvxzuoSiKYgSRkZH83//9H/v27Utat3v3bgBatWpFXFwc\nf/75p17rgoODM/8bMKA0z52klAlCiMHADhKL/nwp5WkhRP/Et+UvQAchxEAgDogB3jNmaEVRrFvO\nnDn59NNP2bRpU9K6/fv3U7lyZQAqVarEH3/8gRBCr3W+vr6Z/00YiF5/AEsptwElX1k3J9nXM4AZ\nho2mKIqW9u/fz+rVq/H390dKSVhYGF9+aVq3gCTvng0PDydbtmwAODs7c+vWLezs7PRaZ85UL6ai\nKKkqUKAAPj4+bNiwwSDHe/jwIcOGDWPevHkpvn/q1Cl27tyJEK93Affo0QMXF5cU99PpdNg+u7qe\nkJCAra2t3uvMmSriiqKkqHbt2nz//ff4+Pjw8OFDnJycDHJcJycnihQp8sb3y5QpQ5ky+o1STl7o\nPTw8iIqKAhJ/UeTJkwcgzXXu7u7p/yZMiCriiqKkKCYmJqlwBwYG0qxZM/bv30+tWrVYvnw5WbJk\noXLlyhw8eJDQ0FBatWrFrl27aN++PatXr6ZZs2acOnWK3r17c/jwYTZu3EjTpk3Zs2cPPj4+nDlz\nhoMHD+Lh4UHTpk2T2n1+Jv4qIQTdu3fH1dU1aV3y7hRfX18OHz5M06ZNCQkJoUGDBtjZ2XHo0KE0\n15kzVcQVRUlRWFgYdevWBRL7jq9cuYK3tzdCCE6dOkXu3LkJDg5mypQpnDt3jqJFi5I3b16yZ89O\niRIlsLGxIWvWrAAULFgQNzc3cubMiaOjI5UrV2bEiBF8/PHHODs7v9SuPmfiUVFRzJ07lzNnzvB/\n//d/9OvXj/r167N161bWrFmDEIJ3330XKSWBgYFprjNnaY4TN2hjapy4RVDjxK3b6tWreeedd7h1\n6xYXL16katWqBAYG4uDgQMeOHTl8+DB+fn4EBgaSN29eWrRowYoVK4iOjkan0xETE0P16tU5d+4c\nZcuWRafTJY0WUVKW2jhxVcSVdFNFXEmP2NhY+vTpw7Rp03Bzc9M6jllKrYir/3aKohiVo6MjS5Ys\n0TqGxVIPhVAURTFjqogriqKYMVXEFUVRzJgq4oqiKGZMFXFFURQzpoq4oihmYeLEiWzcuJHvvvsu\n1e3u37/PqFEvnue+fft2pk+fzowZM4iJiTF2zEyniriiKCYvPXOAL1u2jIiICADu3bvH4sWL+fjj\njwkPD+fMmVef727+VBFXFMXk7d+/n0qVKgEv5gBPyfnz5ylcuHDS8sqVK6levToAY8aMSTqGJVFF\nXFGUFAUHBzN69GgGDBhAly5dNH0CTkpzhackLCwMb2/vpOWTJ09y/fp1AgMDmTp1aqZkzWzqjk1F\nUVLk7u6Os7Mz9evXx8/PD0dHR4McNyQkhGrVqgH6zx2uzxzgBw4coHbt2kRHRyet0+l0uLi4JM2o\nuHXr1pdmTLQEqogripKikiVLcvjwYUaNGoW9vb3BjhsYGJhUxPWdO/zVucJTmgP87NmzXLhwgYiI\nCC5cuMDBgwfJly8f+fPnB8DNzY2TJ0+qIq5YnocP4dChF8tOTlCzpnZ5FNMgpeTp06dJBfzy5ctJ\nc4ePGzeOe/fuceDAAR49ekTjxo2Tvvb19U2aJ9zd3Z2NGzfSrFkz7t+/j5QSIQQPHz4kR44ces8d\nntJc4QBXrlyhUKFCAPTq1Stp3cmTJ6lRowbx8fHs2bMHSLzIWb58eaN/bplNFXGFx4/hzz9fLOfK\npYq4AlevXqVKlSpJy9OmTeOnn37i/PnzZM2alcmTJ/Pzzz8jhODjjz9O+rpXr15J84S7urri5uaG\nq6srZ8+exdbWlp49eyb1b+t7Jp7SXOH379/n/fffZ//+/UnbPXnyhOnTp3Po0CH27dtH3bp12bNn\nDwsWLMDW1pbGjRsb/oPSmJqKVuHGDfjllxfLuXLBxx+/eXs1Fa11mjlzZtLc4Z06deLvv/+mYsWK\n6HQ6jh8/nvT16dOnKVu2LAkJCZw/fz5pHvEePXqwaNEiSpcuTbVq1QzaRWPp1HziSqpUEVeMQc0j\nbjhqPnFFUTKdmkc8c6hx4oqiKGZMnYkrmeb249ucDD/JmTtnOHv3LLejbhOviydBlwCAZw5PirgW\nobBrYarmr0oh10IaJ1YU06eKuGI0UkpC/gth87nNbDm/hUv3L1HeozylcpWiZO6S1PKqhZ2NHXY2\nduikjmsPrnH5/mX2Xd3HoMBB5HDMQcMiDWlWohlNSzTFzkb9uCrKq9T/CsXgop5GsfT4Un4O+Zl4\nXTytS7ZmauOp1PKqhb2tfiMSdFLHyfCT7P53N5P2T2LAlgF8WOlD+lTuQ0GXgkb+DhTFfKgirhjM\no9hHTN4/mVmHZ+Fb0JfpTadTr3C9FG+pTouNsKG8R3nKe5Tn05qfcjL8JHMOz6HSnEq0KdmG8f7j\n8XLxMsJ3oSjmRV3YVDJMRwK/Hp1PyYCSXH5wmZC+IazvvJ76Req/VQFPSdk8ZZnebDoXP7lInmx5\nqDinIiN3jiQyJtIgx1cUc6WKuJIhtzjGXKqy5PgiNnTewJK2Syias6jR2nPN4sr3Db/nxMATRMZE\nUnZWWTae3Wi09hTF1KkirrwVHQkEM4klNKQ6Q/mj2158CvhkWvv5s+dnbqu5LGu3jGHbh9Ht927c\ni7mXae0riqlQRVxJtysPrrAQPy6ynb4cpiI9DNZtkl5+hf04NuAYblncKD+rPAeuHdAkh6JoRRVx\nJV32XdmH74IalKQV3dmNK9qP5c7mkI1pTacxu8Vs2qxow+zDs1HTOyjWQhVxK/D0KURGvng9evR2\nx5l9eDYdV3fk11aLqM1IRCo/Pvfvv9xmQsJbhk+HFu+0YH/v/QSEBNBnYx+exD8xfqOKojG9hhgK\nIZoA/0di0Z8vpfzhDdv5AAeA96SU6wyWUsmQixdh5coXy0WKQI8e+u+foEtg6Lah7L60m+BewRRx\nKcHfaewzdy48m8MfgE8+gcyYA6lErhIc7HOQnut70vS3pmzovIEcjjmM37CiaCTNM3EhhA0QADQG\nvIEuQohSb9huErDd0CEV7cQlxPHB7x9wMuIkB/scpESuElpHSpOzgzMrO6zE290bv4V+3Hqc8vMY\nFcUS6NOdUg04L6W8IqWMA1YArVPY7mNgDRBuwHyKhqLjommzsg1RT6MIfD/QrM5obW1smd50Om1L\ntcX3V18u3ruodSRFMQp9ingB4Fqy5evP1iURQuQH2kgpZwHaDFNQDOpR7COa/tYUt6xurO20lqz2\nWbWOlG5CCMb6jeWzWp/hv8hfFXLFIhnqtvv/A0YlW35jIR8/fnzS1/7+/vj7+xsogmIoT2UMrVa0\n4h23d5jTcg42wryvfw+oOgAbYUODxQ0I6hlEYdfCWkdSlFQFBQURFBSk17b6FPH/gOQzDnk+W5dc\nVWCFSBwsnBtoKoSIk1K+ditd8iKumJ54Ypn/uB2VPQswu8Vssy/gz/Wr0o+nCU+pv6g+e3vuVfOu\nKCbt1RPcCRMmvHFbfYr4IaC4EKIQcBPoDHRJvoGUMuk+ayHEAmBTSgVcMW064llLFxxxYmGbhdja\n2GodyaAGVxtMXEIc9RfXJ7hXMB7OHlpHUpQMS/M0S0qZAAwGdgBhwAop5WkhRH8hRL+UdjFwRiUT\nSCSb6Ec8MfRwXm6xc3d/WvNTupbrSvNlzXn89LHWcRQlw/T6nyql3AaUfGXdnDds29sAuZRMto+J\nhHOCHgRhJxy0jmNU4/zGce3BNTqt7sTGLhst9heWYh0so8NTyZBV5xZxlAV0YRMOZNM6jtEJIZjd\nYjYSyYDNA9Qt+opZU0Xcyu3+dzff/D2C9wnEmbxax8k09rb2rO64mqO3jvLdn99pHUdR3poq4lbs\n/N3zdFnbhVkNVuJO6Zfeu3LlCl999RVBQUHodDqNEhqXs4MzG7tsZNbhWWpOcsVsqSJupR7FPqLN\nyjZ8Xe9rauev99J7x46tw8fHh/DwcIYOHUqhQoUYMWIEd+/e1Sit8eTPnp+1ndby4cYPCQsP0zqO\noqSbKuJWSCd1dF/fndpetelfpX/S+vj4J2zZ8hHr13/Gpk2bmDNnDkePHmXbtm2Eh4fToUMH4uPj\nNUxuHNU9qzPl3Sm0XtFaPVhCMTuqiFuhTQ++5fbj20xvOv2lhzls3z6chw+vMmJEKNWrV09a7+3t\nza+//oqjoyMjR47UIrLRda/QnTal2tB5TWcSdJkwb66iGIgq4lbmPFsJejSHtZ3W4mjnmLQ+MvIS\nYWErad16IU5Orq/tZ2try7Jly9iwYQPLly/LzMiZZlLDScTp4pi4b6LWURRFb6qIW5EHXGMDvRjo\nvoJ82fO99N6+fRPx8fkIJ6dcb9zfzc2NdevWMWzYEG7dOmrsuJnOzsaO5e2XMzd0Ljsu7tA6jqLo\nRRVxK5FAHGvpTA2G8k4W35feu3jxHOfObaJmzU/TPE6FChX43/9+Yv36HkhpeaNW8jrn5bd2v9Fj\nfQ+uP7yudRxFSZMq4lbiD77EERdq83qf9k8/jadGjWFkyfJ6N0pKunb9AFtbB06dWmvomCbBv7A/\nH1f7mPfWvEdcQpzWcRQlVaqIW4G9N7ZwkuW0ZfFrz8U8ceIE+/f/QfXqH+t9PCEE9epNJChoHDoL\nvQg42nc0ORxzMD5ovNZRFCVVqohbuFuPb/HV4T60ZxlO5H7t/QkTJjBw4EgcHJzTddxixRqTNWtO\nTp5cbqioJsVG2LCw9UIWHF1A0OUgreMoyhupIm7BdFJHz/U96Vi0LwXxfe39yMhIdu7cSbdu/VPY\nO3WJZ+PfsHfvBBIstMvBw9mD+a3m0/337mr8uGKyVBG3YAEhAUQ+iWRAma9SfH/z5s3Uq1cPJ6e3\nm/SqSJF65MjhxbFjizMS06Q1LdGUdqXb0XdTXzVRlmKSVBG3UCdun2Divon81u437G3sU9zm999/\np23bthlqp169iezb9zVPnz7N0HFM2aSGk7hw7wLz/5mvdRRFeY0q4hYoNj6Wruu6MrnhZIq7FU9x\nm+joaHbv3k3Lli0z1FbBgrVxdS3Chg3rM3QcU5bFLgvL2i3j892fcynyktZxFOUlqohboHFB4yjm\nVoyeFXu+cZsdO3ZQtWpV3NzcMtyej88gZs+emeHjmDLvPN6Mqj2Knht6orPA8fGK+VJF3MIcvH6Q\nhUcXMrv57JfmRXmVIbpSnitVqg3nzp3l1KlTBjmeqfq0xqfopI5pB6dpHUVRkqgibkGi46Lpsb4H\nAc0CUn0IcHx8HJs3b6ZNmzYGadfW1oHevfswa9YsgxzPVNna2LKw9UK+C/6O0xGntY6jKIAq4hbl\ni91fUCVfFTqU6ZDqdufO7aNYsWJ4enoarO0+ffrx22+/8fixZT98uJhbMSbWm0j39d2J11netLyK\n+VFF3ELsvbyX1adWE9AsIM1tjxwxXFfKc15eXvj5+bFsmWXOcJhc/yr9cc3iypQDU7SOoiiqiFuC\n6LhoPtz4IbOaz8Ita+oXKqXUERq63uBFHGDQoEHMnDnT4sdTCyGY23Iu/zvwP87cOaN1HMXKqSJu\nAcbtGYdPAR9alWyV5rY3b/6Do6MzpUqVMniOBg0aEBUVxV9//WXwY5uawq6FmeA/gd4bequHSCia\nUkXczB367xBLji/h5yY/67X9lSt7KVOmgVGy2NjY0L9/f+bPt46bYgb6DMTOxo6AkLS7sBTFWFQR\nN2NPE57Se2Nvfmr8E+7Z3PXa5+rVYN555/V5VAyla9eurFu3jqdPo43WhqmwETbMbzWfifsmcvHe\nRa3jKFZKFXEz9v2f31PYtTBdynbRa3spJVevBlOihPGKeL58+ahRowYnT/5utDZMSYlcJRhVexQD\ntgyw+GsBimlSRdxMnYo4xfSQ6cxqPivVm3qSu3fvAvb2WcmVy8uo2Xr06MGRI4uM2oYp+bTmp9yJ\nvsOS40u0jqJYIVXEzZBO6ui/uT8T/CfgmUP/sd5Xrwbj5VXbiMkStW7dmuvXD/PQSh5vZmdjx7yW\n8xixcwQRURFax1GsjCriZmhe6DziEuIYUHVAuva7ejWYggWN15XyXNasWSlXrgPHjy81elumokr+\nKnxQ7gOG7RimdRTFyqgibmZuPrrJmD/G8EvLX7C1sU3XvteuZU4RB6hatQfHji3Su584KAhWrXrx\num6GJ/Ff1/ua4KvBbL+wXesoihVRRdzMDNk2hL6V+1Leo3y69ouKCufx49u4u3sbKdnLChWqRUJC\nHDduHNJr+ytX4NSpF69Hj4wc0AiyOWRjVvNZDNwykOg4yx+do5gGVcTNSOD5QEJvhvJV3ZSf1JOa\nq1f34+VVC5t0nr2/LSEEFSp05+hR67nACdCkeBN8Cvjw7b5vtY6iWAlVxM1EdFw0HwV+xIxmM8hq\nnzXd+1+7tj/TulKeq1ChO2FhK0lIsNyn/qRkauOp/BL6C6ciLHtqXsU0qCJuJr7d9y3VClSjcfHG\nb7V/Zl3UTM7VtTC5c5figpX1EefPnp9xfuMYuGWgGjuuGJ1eRVwI0UQIcUYIcU4IMSqF91sJIY4J\nIf4RQoQIIYw/js2KnIo4xZwjc5jaeOpb7R8TE014+Any5/cxcLK0lSvXlRMnfsv0drU2sGpiv/ii\nY9bVnaRkvjSLuBDCBggAGgPeQBchxKuzJ+2SUlaQUlYCPgTmGTyplZJSMmjLIMb5jSN/9vxvdYwT\nJ0Lw8CiP/Vt0w2SUt3cnLlzYysOHDzO9bS3Z2tgyu/lsRu0axd3ou1rHUSyYPmfi1YDzUsorUso4\nYAXQOvkGUsrkl+KdAfUQQgNZcnwJj54+YpDPoLc+RmhoMF5emduV8pyTUy4KFfJjyxbruA0/uSr5\nq/Ce93t8vvtzraMoFsxOj20KANeSLV8nsbC/RAjRBvgecAeaGySdlYuMiWTUrlFs6Lwh1THhf/8N\nO3a8WK5QAVolm5X2xIkQPD17vHWOx49harKenPR285Yr15U1a+bz0Udvn2HDBjh+/MVy48ZQ7bWf\nQtMzsd5ESs8ozcHrB6nhWUPrOIoF0qeI60VKuR5YL4TwBb4BGqW03fjx45O+9vf3x9/f31ARLM6X\nf3xJ65KtqVYg9WolJSQkm9Ja98rfQadOhdK5s35T1b5JQgamzC5ZsiUzZgzg1q1b5M2b962O8er3\naC7XC12yuDC50WQGbRnEob6H0n2DlmKdgoKCCAoK0mtbfYr4f0DBZMuez9alSEoZLIQoKoRwk1Le\ne/X95EVcebMjN46w9vRaTn2UsWFqt2/fJiYmGheXQgZKln729k40bdqaFStWMHToUM1yaKVrua7M\nC53HrMOzGFxtsNZxFDPw6gnuhAkT3ritPn3ih4DiQohCQggHoDOwMfkGQohiyb6uDDikVMAV/eik\njkGBg/iuwXdpPm4tLf/88w+lS1fSe6ZDY+nQoStLl1rPXCrJCSGY2XwmE/ZO4NbjW1rHUSxMmkVc\nSpkADAZ2AGHACinlaSFEfyFEv2ebtRdCnBRChALTgU5GS2wF5oXOw87Gjp4Ve2b4WKGhoZQpUznj\noTKobt363LhxgzNnrPOZlGXcy9CrYi9G7hypdRTFwug1TlxKuU1KWVJKWUJKOenZujlSyl+efT1Z\nSllWSllZSllbSmn5D1k0krvRd/lqz1fMbDYTG5Hxe7Gen4lrzdbWls6dO7Ns2TKto2hmrN9Y9lze\nw59X/tQ6imJB1B2bJuaL3V/wnvd7VMhbwSDHM5UzcYAPPviApUuXWu1djM4Ozkx5dwofBX5EvC5e\n6ziKhVBF3IQcvnGYjec28nW9rw1yvMjISMLDwylUqIRBjpdRlSpVIkuWLPz1l/X+odaxTEfcs7kz\n89BMraMoFkIVcROhkzo+CvyI7xt8j2sWV4Mc8+jRo1SoUAFbW9MY1iaESDobt1ZCCAKaBjBx30Ru\nP76tdRzFAqgibiJ+/edXbIUt3St0N9gx//nnHypV0r4/PLn333+f1atX8/Spdc1smFxp99L0rtib\nkbvURU4l41QRNwH3Yu4x5o8xzGg2wyAXM58LDQ2lcmXT6A9/rnDhwpQuXZrt261rZsNXfeX3FX9c\n+oP9V/drHUUxc6qIm4Av//iSDqU7UCmfYc+aTbGIA1bfpQKJFzn/1+h/DN46mARdBm6HVayeKuIa\nC70ZyrrT65hYf6JBjxsVFcXly5cpU6aMQY9rCB06dGDbtm08ePBA6yiaes/7PVwcXZhzZI7WURQz\npoq4hnRSx+DAwXxT/5sM35n5quPHj1OmTBns7e0NelxDcHNzo379+qxdu1brKJoSQhDQLIDxQeOJ\niIrQOo5iplQR19DS40uJ18XTu1Jvgx/bVLtSnuvWrRtLlizROobmyuYpS9dyXfli9xdaR1HMlMFm\nMVTS58GTB4zeNZr1ndcb9GLmc6Ghofj4vN2TfB4/hs2bXyzHxaW9z9atkHx6ltjY1Ldv0aIF/fv3\n59KlSxQpUuStclqK8f7jKT2jNCH/haQ5Y6WivEqdiWtkwt4JNCvRzGj/aUNDQ996eGFsLBw+/OJ1\n7Fja+xw58vI+8WnckOjg4ECXLl1YvHjxW2W0JC5ZXPi+wfcMDhyMTqrnqSjpo4q4BsLCw1h6fCnf\nN/jeKMePi3vKmTNnKFeunFGObyg9e/Zk0aJF6F6dAN0KdavQDTsbOxb8s0DrKIqZUUU8k0kpGbx1\nMOP8xuGezd0obdy4cZZChQrh5ORklOMbSqVKlXB2diY4OFjrKJqzETYENAtgzB9jiIyJ1DqOYkZU\nEc9kq8JWERkTSf+q/Y3WxuXLx03+LBwSR2f07NmThQsXah3FJFTOV5m2pdoyds9YraMoZkQV8Uz0\n+OljPtv5GQHNArCzMd415StXTphFEQfo2rUrv//+O0+ePNY6ikn4pv43rDq1imO39LgQoSioIp6p\nvt33Lf6F/fEtaNwnz1+9eoLy5csbtQ1D8fDwoE6dOhw5sk7rKCYhl1Muvvb/msFbB1vtlL1K+qgi\nnknO3T3H3NC5TG442ehtmdOZOCRe4AwOXqh1DJPRp3IfouOiWXbCeh+goehPjRPPBFJKPtn6CZ/7\nfk6+7PmM2taTJ/d5/PheqmOvr16FKVNeLGfkSfaG0KJFC3r2HMC9exdwcyuubRgTYGtjS0DTADqs\n7kDLki3J4ZhD60iKCVNn4plgw9kNXH1wlU+qf2L0tsLDT+Ll5Y2NzZv/aRMS4NGjF6/oaKPHSpWD\ngwO1a/fgyJG52gYxITW9atK4WGO+3muYB4QolksVcSOLiYvh0+2fEtAsAHtb489jcvv2CQoVMp+u\nlOf8/fty7NhCEhKsd57xV01qOIlFxxZxKuKU1lEUE6aKuJFNCp5EtQLVqF+kfqa0d/v2cQoWNL8i\nnjfvO7i7e3P69O9aRzEZebLlYWzdsQwOVBc5lTdTRdyI/o38lxmHZvBjox8zrc3wcPM8EweoUqU/\nR9S0rC8Z6DOQuzF3WX1qtdZRFBOlirgRDd02lM9qfYaXi1emtCelJDz8pNkW8dKl2xIREcbdu+e0\njmIy7GzsCGgawPAdw3kU+0jrOIoJUkXcSDad3cS5u+cYVnNYprX58OE17O2z4uJinNv5jc3W1oEK\nFXpy5MgvWkcxKXUK1aF+kfpM3GfYB4colkEVcSOIiYthyLYhTG86HQdbh0xr9/btE+TJYx5n4eHh\n8N9/L17Pp66tUqUvx44tJjb2iUHbi4p6ub179wx6eKOb3HAyC44u4HTEaa2jKCZGjRM3gh/2/0DV\n/FVpVKxRprYbHm4+Rfz33+HmzdfXu7kVJ2/eCmzfvoaKFT8wWHtnzsCmTS+Wvb2hY0eDHd7oPJw9\n+KruVwzeOphd3XYhkk/erlg1dSZuYBfvXSQgJICfGv+U6W3fvn0cDw/zKOKpqVbtE5YunaZGZLxi\nkM8g7kTfYVXYKq2jKCZEFXEDklLyybZPGFl7JJ45PDO9fXM6E0/NO+805+HD++zfv1/rKCbFzsaO\nGc1mqIucyktUETegDWc3cCnyEkNrDM30thMSnnLv3gXc3U3v6fbpJYQN3boNYerUqVpHMTm+BX1p\nVKwR44PGax1FMRGqiBtI1NMohmwbwszmMzP1YuZzd+6cxcWlIPb2WTO9bWNo06Yne/fu5dKlS1pH\nMTmTG05myfElnLh9QusoiglQRdxAJu6bSN1CdfEv7K9J+5bSlfJctmzO9O7dm59//lnrKCbHPZs7\nX9f7mkGBg9R1A0UVcUM4FXGK+f/M53+N/qdZhvDwkxZVxAEGDx7MokWLePjwodZRTE7fyn2JjY9l\n8TH1oGlb6bZEAAAgAElEQVRrp4p4BkkpGbRlEOP8xpHXOa9mOcLDT1jEyJTkChYsSKNGjZg/f77W\nUUyOrY0ts5rPYtSuUdyLMbNB74pBqSKeQUuOL+HR00cMrDpQ0xyJZ+JlNc1gDMOHD2fq1KnExanZ\nDV9VJX8VOnl3YvSu0VpHUTSkVxEXQjQRQpwRQpwTQoxK4f33hRDHnr2ChRCWdUr4Bvdi7jFy50hm\nN5+NrY2tZjmioh7x+PFtcuYsplkGY6lWrRqlS5dmz55FWkcxSRPrTWTL+S0cuHZA6yiKRtIs4kII\nGyAAaAx4A12EEKVe2exfoK6UsgLwDWAVs/uP2jmKjmU64lPAR9Mcly6F4e5eGhsNf5EY09ixY1m7\n9jsSEuK0jmJyXLK4MLXxVPpv7k+c+nyskj5n4tWA81LKK1LKOGAF0Dr5BlLKg1LKB88WDwIFDBvT\n9ARfDWbrha18U/8braPw77+W2ZXyXO3atfHwKMrx40u1jmKSOpbpSIHsBfi/g/+ndRRFA/oU8QLA\ntWTL10m9SPcBtmYklKl7mvCUAZsHMLXxVFyyuGgdh4sXT+LubrlFHKBTp7H8+ee36HTxWkcxOUII\nZjSbwQ/7f+Dy/ctax1EymUEnwBJC1AN6Ab5v2mb8+PFJX/v7++Pv72/ICJnip79+oqBLQTqU6aB1\nFCDxTLxkySZaxzCqsmX9yJGjACdOLKdChW5axzE5xdyKMazmMAZtGcSW97eoCbLMXFBQEEFBQXpt\nq08R/w8omGzZ89m6lwghygO/AE2klJFvOljyIm6OLty7wI8HfuRQ30Mm8x/l4sUT1Klj2WfiAH5+\n49iyZSDlyr0PWGb/f0Z8Vuszlp9czqqwVbxX9j2t4ygZ8OoJ7oQJE964rT5F/BBQXAhRCLgJdAa6\nJN9ACFEQWAt0k1JeTH9k8yClZMDmAXzu+zlFchbROg4A4eHhxMXFkj27aV+G+O03sEv203b3bvqP\nUbhwPZyd83L06AKaN++ToTwXLsCsWS+WvbygRYsMHfI1yY8P0Ls3ODoato3kHGwd+KXFL7Rf1Z53\ni71Lzqw5jdeYYjLS7BOXUiYAg4EdQBiwQkp5WgjRXwjR79lmXwFuwEwhxD9CiBCjJdbQkuNLuBdz\njyE1hmgdJUlYWBjFipUzmb8K3uTuXbh9+8Ur/i26toUQNGr0I3v2jCU6+nGG8sTGvpzn/v0MHS5F\n4eEvt6HTGb6NV9X0qknbUm0Zteu1kcCKhdJrnLiUcpuUsqSUsoSUctKzdXOklL88+7qvlDKXlLKy\nlLKSlLKaMUNrISIqghE7RzC35VzsbEznWRonTpygaFHL70p5rkABH4oUqc/SpZO1jmKyvmvwHYHn\nA9l3ZZ/WUZRMoO7Y1NOwHcPoVr4bVfJX0TrKS06ePGlVRRygQYPvWLNmBtevX9c6iklyyeLC9KbT\n6bupL0/iDfuYO8X0qCKuh63ntxJ8NZgJ/m++uKCVkydPUqyYVdwgm8TFpSBt2/ZnzJgxWkcxWW1L\nt6VcnnJ8vfdrraMoRqaKeBoexj5kwJYBzG05l2wO2bSO8xIppVWeiQN07z6aHTt2cOTIEa2jmKyA\nZgHMC53HPzf/0TqKYkSqiKdh9K7RNCzSkIZFG2od5TVXr14le/bsuLi4aR0l02XLloOJEycyaNAg\nEhIStI5jkvI652Vyo8l8uPFD4tVNUhZLFfFU7Luyjw1nNzCl8RSto7wkJgYiI+HgwZOULFmWmBit\nExleVFTi9/j8FRv7+ja9e/fGycnJIA+OiIt7ub1HRniE5YMHL7eRGb97elToQW6n3Ew5YFo/w4rh\nmM4wCxMTExdDn419mNFsBq5ZXLWO85KQENizB4KDTxITU5a9e7VOZHibN6e9jY2NDXPnzqVGjRq0\nbt2aokWLvnV7V67AtGkvlj09oU/GhqK/Zvbsl5eHDgVXI/9oCSGY02IOPnN9aF2qNaVyvzp3nWLu\n1Jn4G4zdM5ZK+SrRplQbraO8kSU+CCK9ihcvzsiRI+nXr596VNkbFMlZhAn+E+i1oRcJOtX1ZGlU\nEU/BgWsHWHpiKQFNA7SOkqrE52pa30XNVw0bNozIyEgWLFigdRSTNdBnIFnssjD14FStoygGpor4\nK6Ljoum5vicBTQNwz+audZw3Skh4yt2753B399Y6iubs7OyYP38+o0eP5t9//9U6jkmyETbMbzWf\nScGTOHPnjNZxFANSRfwVY3aPoWr+qrQv017rKKmKiDiNq2sR7O2zah3FJFSsWJExY8bQqVMnYlO6\nCqpQNGdR1a1igVQRT2bflX2sDFvJ9KbTtY6Sptu3j5E3bwWtY5iUTz75hEKFCjF8+HCto5isgT4D\nyWqXlR8P/Kh1FMVAVBF/5lHsI3pt6MWs5rPI5ZRL6zhpunXrGB4eqognJ4Tg119/Zdu2baxcuVLr\nOCbJRtiwoPUCfvzrR47fPq51HMUA1BDDZ4ZtH4ZfIT9al2qd9sYm4PbtY9Ssqc44n7t2Df77D8CF\nr75axUcfNaF8+fKULl3aKO09fQqhoS+WbW3BR9tHrXLpUuJsic95eUGBFGYoLuRaiP81+h8frPuA\nQ30P4WhnxPlxFaNTRRzYeHYjuy7t4tiAY1pH0YuUUnWnvOL8ediXNGlfZXr3nkLTpk3Zv38/BVKq\nZBn05Als2/Zi2cFB+yJ+6hQcOvRiuUGDlIs4JN4EtOHsBsbuGcsPjX7InICKUVh9d0p4VDj9N/dn\ncZvF5HDMoXUcvdy9exMAZ+d8GicxXQ0bdmPQoEE0adKEyMg3PmjKaj2/CWjx8cUEXw3WOo6SAVZd\nxKWU9NvUjx4VelCnUB2t4+jt4sVjeHiUN/kHQWhtxIgRNGzYkNatW/PkiQXOTZBBebLlYU6LOXT7\nvRsPnjzQOo7ylqy6iM8Lncfl+5dNcorZ1Fy4oC5q6kMIwZQpU/D09GTkyPbExUVrHcnktCrZisbF\nGvNR4EdaR1HektUW8dMRp/nijy9Y3n652V3YSTwTV0VcHzY2NixatAhX19wsWdKImJh7WkcyOT81\n/okjN4/w2/HftI6ivAWrLOKx8bF0WduFb+t/S2l344xeMKaLF9VFzfSwt7dn/PiFeHrWZMGCOjx8\nqJ4IlJyTvRPL2i1j6PahXIq8pHUcJZ2ssoiP3jWaYm7F6Fu5r9ZR0i0mJoabNy+RO7f5/fLRko2N\nDe+++yMVKvTk119rc/36Qa0jmZRK+Srxue/ndF3XlbiEOK3jKOlgdUMMt57fytrTazk64KhZXhgM\nCwvD07MEdql0AV24AIsWvViOisqEYGaidu0R5MpVghUrWlO9+lB8fUchxOvnMhERL3+GLi7QxsgT\nWm7cmDjPeGqioiKJiXmITpdAlSo6SpbMSv78+Q3yszy0xlB2XNzB2D1j+b7h9xk+npI5rKqIX394\nnV4berGq4yrcsprn03COHTtGsWKpd6U8fpz4UlJWqlQb8uWrwrp1Xbl0aTetWy/AxcXrpW1iYxNv\nnnnOPRPmQvvvvxc36+h08dy+fZyrV/fz339/c+/eee7ePY9OF0/WrDkRwobs2W2JjX3EkydP8Pb2\nxsWlHLlzd6Bo0QYp/mJKi42wYXHbxVSeUxn/wv40Lt7YwN+hYgxWU8TjEuLovKYzQ6oPoW6hulrH\neWvHjx9Ps4graXNx8aJHjz8IDp7EnDkV8fH5iNq1R+Lg4KxJHiklN2+eISRkOxcv7uDq1WBcXLzw\n8qpNkSL18fH5CDe34jg55U46627XDsqXhzt37hAWFsaSJUfYvHkEsbEPqVTpQypX7gvkSVeOPNny\n8Fu73+i8tjNH+h0hf/b8RvhuFUOymiL+1Z6vyO6YnVG+o7SOkiHHjh2jZcsWRnl8mLWxsbGjbt0v\nqVChO7t3f0FAQEnq1h1LhQrdsLd3Mnr7jx49IihoN1u3bmXbtm08fgyFCzemUqUPadduKVn1/Gsx\nd+7c+Pn58fixHwUKfMrNm6EcOTKHHj3KMnXqZHr06JGu7ha/wn4MqjqI99e+z67uu7CzsZoyYZas\n4sLmlnNbWHZiGUvaLsHmLf7MNBVSSr26U5T0cXEpSLt2S3nvvfWcP7+ZqVO92Lp1CBERpwzajk4X\nz/XrB9m37xsWLvSnTJn8zJgxg5IlS7Jt2zYmTrxMy5a/UKZMe70L+KuEEOTPX4WWLX/hxx+3M23a\nNJo2bcqVK1fSdZwv6nyBva0944PGv1UOJfNY/K/Yy/cv8+HGD1nTaQ25nXJrHSdDrl69iqOjIzlz\npu9PZEU/BQr40KXLJu7fv0Jo6DwWL26Io2N2ihZtRKVKjQgPr4m7u7teZ7VxcXHcufMvERGnuX79\nb65f/4ubN4/g6lqEokUbUbv2KCZPrkuBAtmS9gkKMuz38847lQgJCeHHH3/Ex8eHNWvWULeufl2J\ntja2/NbuN6r+UpXqBarTsmRLw4ZTDMaii3hMXAztVrZjtO9ofAv6ah0nww4fPkzVqlW1jmHxXF0L\nUb/+ROrVm8CtW8f499+d7N0bQKlSvdDpdBQvXhxPzyJcvpwNW1tH7Owc0eliOHLkHpGRkdy4cYMr\nV66QNWsBcucuRf78Pvj6fo6nZ3WyJHvodrZsqYQwEHt7ez7//HN8fHxo3749CxcupHnz5nrtmydb\nHlZ1XEWr5a048OEBirsVN3Ja5W1YbBGXUjJwy0BK5i7JkOpDtI5jEIcOHaJatWpax7AaQtiQL18l\n8uWrRK5cI+nXD+7du8eFC+c5c+YyW7bEEB//hPj4WLJkyUqXLm64ubnh4eFB8eLFmTTJkdSe3fz0\naeIomOfS+5zn+PiX909I5WE9DRs2ZPPmzbRq1YqpU6fy/vvv69VGDc8ajPcfT/tV7fnrw79wyoRr\nBUr6WGwRn314NkduHuHghwfNcjx4SkJCQhgxYoTWMazS3bvw/fcAbkB1oDoVK75438EBOnRI3zFn\nzsxYpo0bE1/6ql69On/88QdNmjQB0LuQD6w6kL+u/5U026el/H+yFOZ7lS8VB64dYFzQONZ1Wkc2\nh0z4mzUT6HQ6jhw5go/Wk1YrZs3b25vAwECGDBlCSEiIXvs8n7Y2LDyMqQenGjmhkl4WV8SvPrhK\nh1UdWNhmISVyldA6jsGcO3eOXLlykTu3eV+cVbRXrlw55s+fT7t27bh+Xb95ZJzsnVjfeT0/HviR\nbRe2pb2DkmksqohHPY2i1fJWDK85nGYlmmkdx6BCQkLUWbhiMK1atWLw4MG0adOG6Gj9pugt6FKQ\nVR1X0f337py9c9bICRV9WUwR10kdPdb3oFK+SgyrOUzrOAZ36NAhVcQVgxo1ahSlS5dmwIABeu/j\nW9CX7xp8R6sVrbj/5L4R0yn6spgiPj5oPDcf32R289kWeeFFjUxRDE0IwZw5c/j7779Zs2aN3vv1\nqdyHpsWb0n5Ve54mPDViQkUfehVxIUQTIcQZIcQ5IcRr960LIUoKIQ4IIZ4IITL9NHjBPwtYenwp\n6zqtM7sHPOjj6dOnnDhxgsqVK2sdRbEwTk5OLF68mMGDB3Pr1i2995vy7hSy2WdjwOYByPSOjVQM\nKs0hhiJxOrQAoAFwAzgkhNggpTyTbLO7wMeAkSfrfN2uf3cxevdo9vbci4ezR2Y3nylOnDhB0aJF\ncXbWZnImU3TjBhw//mI5PPzl9yMiXn5fz+t3Fi08HJLX6Vy5oECBxKGHffr0oVevfkyatOGNf8lm\nywbFiiV+bWtjy7L2y/Bb6Mf3wd/zRZ0vMuE7UFKizzjxasB5KeUVACHECqA1kFTEpZR3gDtCiBZG\nSfkGJ26f4P2177Om0xpK5S6VmU1nKtUf/rpjxxJfb3L+fOJLeeHcOdi168Wyj09iEQcYO3Ys5ctX\nZ+zYBVSq1DvF/QsWfFHEAZwdnNnUZRM159ekiGsRupTrYsT0ypvo051SALiWbPn6s3WauvbgGs2X\nNWdak2lmPbWsPlQRV4zNwcGBb79dzK5do9L1+Lr82fOzuctmhm4fys6LO42YUHkTs7yweSf6Du8u\nfZch1YdYxW//kJAQdVFTMboSJcpRpcoAduwYnq79ynmUY03HNXRd15VD/x0yUjrlTfTpTvkPKJhs\n2fPZurcyfvz4pK/9/f3x9/dP1/6PYh/R7LdmtCnZhuG10vfDZo6ioqL4999/KVeunNZRFCtQp87n\nzJzpzb//7qJo0Yb671eoDvNazaPVilYE9QiiZO6SRkxp+YKCggjSc1pLfYr4IaC4EKIQcBPoDKR2\n+pvq+L7kRTy9YuNjabuyLRU8KvBdg+/e+jjmJDQ0lLJly+Lg4KB1FMUK2Ns70aTJNAIDBzNw4HFs\nbfX/uWtVshV3ou/QeGlj9vXaR0GXgmnvpKTo1RPcCRMmvHHbNLtTpJQJwGBgBxAGrJBSnhZC9BdC\n9AMQQngIIa4BnwJjhBBXhRAGHUoRlxDHe2vewyWLC7NazLLIseApUV0pSmYrWbIVuXKV4K+/fkr3\nvr0r9WZI9SE0WNyAG49uGCGd8iq9ZjGUUm4DSr6ybk6yr28DXq/uZyjxunjeX/c+CTKBVe1XWdXj\novbt20eXLpbf76+YliZNpjF3bjXKlXsfl3SeUX9a81Ni4mNouLghQT2DyJNNPcTEmEy+GiboEuj+\ne3ceP33M+vfW45COP+/MXUJCAvv27WPIkF9YvfrF+lfHRCvK27hwgZd+ru4nu4s+Z86iVKs2mF27\nRtG+/fJ0H/uLOl8QHRdNoyWN2NNjD25v+bg5JW0mXcQTdAn03tibiOgINnbeaJF3Y6bm6NGj5M+f\nn7g4D8LCtE6jWJrIyMTXm9SqNYKAgJJcu/YXXl410338ifUmEhsfS/1F9dnZbSfu2dwzkFZ5E5Md\nYhiXEEfXdV258egGGzpvIKt9Vq0jZbo//viDevXqaR1DsVIODtmoX/9bduwY9la31gshmNxoMs1L\nNKfeonrceqz/bf2K/kyyiMfGx9JxdUei4qLY1GWT1T4Sas+ePdSvX1/rGIoVq1ChGwkJTwkLW/VW\n+wsh+Kb+N3Ty7oT/Qn/+e/jWo5OVNzC5Ih4dF03rFa2xt7Vnbae1ZLHLonUkTcTFxREcHIyfn5/W\nURQrJoQN7747hV27RhEX9+QtjyEY6zeWnhV7UndhXS7eu2jglNbNpIr43ei7NFjcAA9nD5a3X25V\nFzFfdfjwYYoVK0auXLm0jqJYucKF/cmbtyI7d07L0HFG+45mRK0R1F1Yl39u/mOgdIrJFPEr969Q\n+9fa+BXyY2HrhVY1jDAle/bsUf3hislo1Ggy27f/j4iIiAwdZ0DVAUxrMo3GSxsTdDnIMOGsnEkU\n8eO3j+O7wJeBVQcyqeEkq7mRJzXqoqZiSnLleofq1d/n66+/zvCxOpTpwIoOK+i0uhPLT6R/+KLy\nMs2L+OZzm2m4uCE/NvqRITWGaB3HJMTGxvL3339Tt65lz86omJeWLceyfPlyzp7N+PM16xepz67u\nic8C+Hrv1+rBEhmgWRGXUvLTXz/Rb1M/NnbZyHtl39Mqisk5ePAgpUuXxsXFResoipIke/bcjBw5\nktGjRxvkeOU9ynPww4NsPreZ7uu7Exsfa5DjWhtNinhsfCz9NvVj0bFFHOxzkBqeNbSIYbLU0ELF\nVH3yySeEhoby559/GuR4+bLnI6hnEDFxMfgvUkMQ30amF/FrD65RZ0Ed7sbcJbhXsJrpLAXqoqZi\nqrJkycJ3333H8OHD0el0Bjmmk70TqzquokWJFvjM9eHPK4b5BWEtMr2IV5tXjY5lOrK201qyO2bP\n7OZN3t27dzl69Ci+vr5aR1GUFHXp0gUpJcuXG+6ipI2wYUzdMcxvNZ8Oqzvw898/q35yPWV6Ef+t\n3W+MqD1CjUB5g99//53GjRuTLVs2raMoSopsbGyYNm0ao0ePJioqyqDHblqiKQd6H2DB0QV0WN2B\nyJhUJndRAA2KeP0iqq83NatWraJTp05ax1CUVNWqVYs6derwww8/GPzYxdyKcfDDg3hm96TSnEoc\nuHbA4G1YEuu+o8ZE3LoFjx/D3bsRHDwYwpQp67lw4cX7yacIVUyTTsdL/2YAWvcG3L37cqa7dw17\n/B9++IGKFSvSu3dvChcubNBjO9o5Mq3pNBoUbUDblW0ZWHUgY+qMwd7W3qDtWAKRmf1OQgip+rle\nt2oVnDoFhw/P5sqVvW81f7OiGFvBgtC798vrvv76a06ePMmqVW83QZY+/nv4H3039eV21G0Wt1mM\ndx5vo7VlqoQQSClT7IPW/GYf5YWwsFWUKaO6UhTz8dlnn/H333/r/VDft1EgRwG2vL+FgVUH4r/I\nn0nBk4hLiDNae+ZGFXET8fjxLW7eDKVEiaZaR1EUvTk5OTFlyhQGDRpEbKzxbtYRQtCnch8O9T1E\n0OUgqvxShYPXDxqtPXOiiriJOHVqLe+80wI7K516VzFf7du3p3jx4kyaNMnobRV2LczWrlv5os4X\ntFvZjkFbBnEv5p7R2zVlqoibiLCwlXh7q6kHFPMjhGDmzJkEBARw6tSpTGmvc9nOhA1KfGZhqYBS\nzAiZQbwu3uhtmyJVxE3AvXv/ER5+kmLF3tU6iqK8FU9PT8aPH0/fvn0NdidnWnJmzcnM5jPZ3X03\n686so+Lsimw9v9XqbhJSRdwE7Nw5G2/vTthZ2YOgFcsycOBApJTMmTMnU9st51GOXd128U39bxi+\nYzh+C/3Yf3V/pmbQkhpiqLH79+9TsGBxevUKIWfOolrHUZQ3cnaG0qVfLGfPDslnS46OhkWLTjF6\ntB8//vgX+fMXf+0YzZu/vLxtGyQkvFiuXx+yZuCZ6PG6eJYeX8r4oPF45/HmyzpfUtOr5tsf0ESk\nNsRQ3eyjsRkzZlCxYjNVwBWT9/gxHDr0YtnD4+Ui/vQp3L5dhtq1xzN2bCc+/PDASxfqhXi9iIeG\nJu73XJ06GSvidjZ29KzYky5lu/DrP7/SZW0XiuYsypg6Y6hfpL5FTvehulM09PjxY6ZNm0abNp9r\nHUVRDMbHZxBubsXYvn24Zhkc7RwZ6DOQ8x+fp3uF7nwU+BFVfqnC4mOLLW7eclXENTRnzhz8/f3x\n9Cyd9saKYiaEELRsOY+LF7cRFma8Ozn1YW9rT8+KPTn10Sm+qf8NS44vofC0wowPGs+1B9c0zWYo\nqohrJCYmhilTpjBmzBitoyhGFBv7iGvX/iI29pHWUTJVliwudOiwisDAwdy9e07rONgIG5qVaMbO\nbjvZ2W0nd6LvUHFORVoub8nGsxvN+g5QVcQ1Mn/+fKpWrUqFChW0jqIYSWzsIxYsqMPChXVZsKCO\n1RXy/Pmr0LDhDyxZ8i4PTOist2yesgQ0C+Dq0Ku0L92eyfsnU+CnAgwOHMzB6wfNbohipl/YtMQL\nCxmhPg/rEBFxioiIMDyt7FGElSr1IibmHkuWNKJXr31AHq0jJcnmkI2eFXvSs2JP/o38l2UnltFz\nfU+exD+hXel2tC/dnppeNbERpn2um+lF3Nx+yxnakydPqF27Nj169OCTTz4BXsxiqFiW52fit28f\nw929DO7u1jf7HkCtWsOJjX3A0qWNGTNmD66urlpHek3RnEX5su6XjKkzhpPhJ1l7ei0DtgwgIiqC\nZiWa0bxEcxoVa0QOxxxaR32NGieeCZYtg5s3E79et+4jHj26Tffuq5POwmNiIN467xi2eLGxj5g0\nKQejRz/E0cIeR2hjA8kfQKXTwZse9COlZNu2ody4cYCePX/H1dUTgEev9DANGwY5TKhOXrx3kS3n\ntxB4PpD91/ZTMW9FGhRpQIMiDajuWR0HW4dMyZHaOHFVxDPBvHlw/TqcPLmSP/74gn79QsmSxUXr\nWEommTBBMG6c9f3cv0pKyf79P/D33z/TocMKChWq+9o2plbEk4uOiyb4ajC7/93N7ku7OXPnDFXy\nV8HXy5faBWtTrUA1cjvlNkrb6mYfE3Du3Ba2bh3MBx9sVwVcsUpCCHx9R5M3byVWr+5InTpjqFbt\nY7O5LuRk78S7xd7l3WdzHD2MfcjB6wcJvhrMlL+mcPjGYXJlzUW1AtWonK8yFTwqUDFvRTycPYya\nS68zcSFEE+D/SBzNMl9K+dqD9YQQPwNNgSigp5TyaArbWN2ZeEJCAm3aTGDv3l/p0GEFBQuqp9hb\nG3Um/rrIyH9Zs+Y9hLDl3Xd/TPp/Ycpn4mnRSR3n7p4j5L8Qjt46mvSyt7WnjHsZyuQuQ2n30pTM\nVZISuUrglcMLWxtbvY6doe4UIYQNcA5oANwADgGdpZRnkm3TFBgspWwuhKgOTJNSvnYZ3lSKeFBQ\nEP7+/kZv5+bNm/Tq1YtLl57QsuUKnJ3zGr3N9Lp8OYjChf21jmESjPVZmGMRz4yfCyl1HD/+G3v2\nfEm+fFWoW/crJk+uiIuLaZ2ZZ6ReSCm58egGpyJOJb3O3zvPubvnuBtzl0IuhSiSswhFXItQ2LUw\nXjm88HLxwiuHF/my50vqc89od0o14LyU8sqzg60AWgNnkm3TGlj8LPTfQggXIYSHlPL2W33nRmbM\nIp7Y77efmTNnsnXrVgYOHEjbtl9z86Zp9lypIv6C+ixeyIzPQggbKlToRpkyHQgJmc7KlW3Zty8b\n3bp15b333qNo0aIm0dWSkXohhKBAjgIUyFGARsUavfRe1NMoLt+/zKX7l7gUeYkrD65w5OYRrj24\nxrWH17j9+DY5HHOQP3v+VNvQp7IUAJKP1L9OYmFPbZv/nq0zySJuKHFxcURERHDlyhWOHDnC4cOH\n+euvvwAYNGgQM2fOxNXVlXnzNA6qKCbM3j4rtWuPpFatz6hV6wDr1/+Gr68vOp2OatWq4ePjwzvv\nvEPBggXx8vIib9682Nub/1PvszlkwzuP9xsf/KyTOu5E3+Hmo5tUHFTxjcfJ9NPDli1bZnaTrzl7\n9ixHjhwBXh63LqVMWpZSotPpSEhIICEhgfj4eJ48eUJMTAwxMTHcv3+fhw8fkjt3bjw9PalUqRK1\narm0d1kAAANoSURBVNXik08+oWLFitjYvLhBIFeul6fbNCXOzpAvn9YpTIMxPwtz+4y1+bmwoXZt\nXxo39mXmzJlcv36dkJAQDh06xNq1a7l27RpXr14lPDwcR0dHcuTIQY4cOciSJQsODg44OjpiZ2eH\njY0Ntra22NjYIIR46ZVc8uXUzviT1wtTpE+feA1gvJSyybPl0YBMfnFTCDEb2COlXPls+Qzg92p3\nihDCvDoGFUVRTERG+sQPAcWFEIWAm0BnoMsr22wEPgJWPiv691PqD39TCEVRFOXtpFnEpZQJQojB\nwA5eDDE8LYTon/i2/EVKGSiEaCaEuEDiEMNexo2tKIqiQCbfsakoiqIYlmlPz2VkQojhQgidEMJN\n6yxaEUJMFkKcFkIcFUKsFUKY6a0Wb08I0UQIcUYIcU4IMUrrPFoRQngKIf4QQoQJIU4IIT7ROpPW\nhBA2QohQIcRGrbO8idUWcSGEJ9AIuKJ1Fo3tALyllBWB84BVPSvu2c1sAUBjwBvoIoQopW0qzcQD\nw6SU3kBN4CMr/iyeGwKY9ByjVlvEganACK1DaE1KuUtKqXu2eBDw1DKPBpJuZpNSxgHPb2azOlLK\nW8+ny5BSPgZOk3i/h1V6dqLXDDDpOz2ssogLIVoB16SUJ7TOYmJ6A1u1DpHJUrqZzWoL13NCiMJA\nReBvbZNo6vmJnklfODTNe8ENQAixE0g+fZgg8R/jS+ALErtSkr9nsVL5LMZIKTc922YMECelXKZB\nRMWECCGcgTXAkGdn5FZHCNEcuC2lPCqE8MeEa4TFFnEpZaOU1gshygKFgWMi8TYtT+CIEKKalDI8\nEyNmmjd9Fs8JIXqS+Gdj/UwJZFr+AwomW/Z8ts4qCSHsSCzgS6SUG7TOo6HaQCshRDMgK5BdCLFY\nStld41yvsfohhkKIS0BlKWWk1lm08Gya4SlAXSnlXa3zZDYhhC1wlsRZOm8CIUAXKeVpTYNpRAix\nGLgjpRymdRZTIYTwA4ZLKVtpnSUlVtkn/gqJCf+plAmmA87AzmdDqWZqHSgzSSkTgOc3s4UBK6y4\ngNcGugL1hRD/PPt5aKJ1LiV1Vn8mriiKYs7UmbiiKIoZU0VcURTFjKkiriiKYsZUEVcURTFjqogr\niqKYMVXEFUVRzJgq4oqiKGZMFXFFURQz9v+aGmW6VhQRgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7328ceee50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# test\n",
    "mu = -1.0\n",
    "sigma = 2.0\n",
    "nsamp = 1000\n",
    "xsample, ysample, accept, c, normalize = sampling(mu, sigma, nsamples=nsamp)\n",
    "plot_samples(mu, sigma, xsample, ysample, accept, c, normalize, trace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n"
     ]
    }
   ],
   "source": [
    "# animation plot\n",
    "# make an animation\n",
    "mu = -1.0\n",
    "sigma = 2.0\n",
    "nsample = 1000\n",
    "xsample, ysample, accept, c, normalize = sampling(mu, sigma, nsamples=nsample)\n",
    "divisor = 10\n",
    "\n",
    "for i in range(nsample):\n",
    "    plot_samples(mu, sigma,\n",
    "                xsample[0:i], \n",
    "                ysample[0:i], \n",
    "                accept[0:i], \n",
    "                c, 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": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkz\nODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2Nj\nY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQED\nEQH/xAAbAAEAAwADAQAAAAAAAAAAAAAAAQQFAgMGB//EAEwQAAEDAQIHCQsKBgIDAQEAAAABAgME\nBRESFSExVXHRBhM1QVFyk5SxFCIyMzQ2YYGRweEWI0JSVGRzobLCJFNigpLidINDY/BEJf/EABkB\nAQEBAQEBAAAAAAAAAAAAAAABAgMEBf/EACcRAQEAAQMCBgIDAQAAAAAAAAABAgMRMQQhEjIzQVFx\nFIEiNGET/9oADAMBAAIRAxEAPwD5+AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAD7hiKx9FUPV2bBiKx9FUPV2bC+VKq0aekkwJVerkbhuw\nGK7BbyrdmQDrxFY+iqHq7NgxFY+iqHq7Nhbknjjp3Tud821uGqplyXXnXS1kdUirGyVqXIt741bf\n7QK+IrH0VQ9XZsJxFY+iqHq7NhWdDJWW5VxLV1EUcUUStbE/BS9cK/sQsYq+/wBd0wE4isfRVD1d\nmwYisfRVD1dmwjFX3+u6YYq+/wBd0wE4isfRVD1dmwjEVj6Koers2DFX3+u6Y6qmzHsjRWV9cq4T\nU8bxKqX/AJAd2IrH0VQ9XZsGIrH0VQ9XZsIxV9/rumGKvv8AXdMBOIrH0VQ9XZsGIrH0VQ9XZsIx\nV9/rumGKvv8AXdMBOIrH0VQ9XZsGIrH0VQ9XZsIxV9/rumGKvv8AXdMBOIrH0VQ9XZsGIrH0VQ9X\nZsIxV9/rumGKvv8AXdMBOIrH0VQ9XZsIxFY+iqHq7NgxV9/rumMaujqILTbAy0KzAVW/+XlN4YXO\n7QbOIrH0VQ9XZsJxFY+iqHq7Nh0T2a9joUZXVy4UiI753iuU7sVff67pjAnEVj6Koers2DEVj6Ko\ners2HHFX3+u6YnFX3+u6YCcRWPoqh6uzYMRWPoqh6uzYRir7/XdMMVff67pgJxFY+iqHq7NgxFY+\niqHq7NhGKvv9d0wxV9/rumAnEVj6Koers2DEVj6Koers2EYq+/13TDFX3+u6YCcRWPoqh6uzYRiK\nx9FUPV2bBir7/XdMdEtmyNqYWNrq7AdhYS776MgFjEVj6Koers2DEVj6Koers2EYq+/13TDFX3+u\n6YCcRWPoqh6uzYMRWPoqh6uzYRir7/XdMMVff67pgJxFY+iqHq7NgxFY+iqHq7NhGKvv9d0wxV9/\nrumAnEVj6Koers2DEVj6Koers2EYq+/13TDFX3+u6YCcRWPoqh6uzYMRWPoqh6uzYRir7/XdMMVf\nf67pgGIrH0VQ9XZsJxFY+iqHq7Nh0LZr+7Ej7vrsDe779947zuxUukK7pgJxFY+iqHq7NgxFY+iq\nHq7NhGKl0hXdMMVLpCu6YCcRWPoqh6uzYMRWPoqh6uzYRipdIV3TDFS6QrumAnEVj6Koers2DEVj\n6Koers2EYqXSFd0wxUukK7pgJxFY+iqHq7NgxFY+iqHq7NhGKl0hXdMMVLpCu6YCcRWPoqh6uzYR\niKx9FUPV2bBipdIV3THS2zHrUyMWvrsBGtVPnePLeB34isfRVD1dmwYisfRVD1dmwjFK6QrumGKV\n0hXdMBOIrH0VQ9XZsGIrH0VQ9XZsIxSukK7phildIV3TAMRWPoqh6uzYMRWPoqh6uzYVLTo5KOhf\nPFX1ivYrbsKW9M6G2AMaqlSjtOsfLG9WzwNSNWsV2EqYXe5OPKntNogDMpWMhsaOhq8PCjpWpLc1\nVyXXLcqalOuzXXV7o6SSaSjSJL98vVGvvyXKuXMa9wuAyImzLuhr95exvzMN+ExXfX9KGu2+5MLP\nx3GbS+cNofgw/vNMAAABBIAgkAAAAAAAAAAebtThxuth6Q83anDjdbD09N5r9JXoyvaDnMoJ3NVU\ncjFuVCwVrS4On5inDDzRVKwJpJo5lle59ypdeprGLub8VPrQ2jprzbUuyTgJIJOKgAAAAAQSAIJA\nAAAAAAAAAAADrlljhbhSvaxFW69VuJjlZKirG5HIi3LdxKdNZTd0sYy9ERJGudfxoi33fkVVZVUb\nJnRI2R80qvyIubbciJ6gNIGclTWuV6JFcuS69i8bruy5Ti6rq2J37blVWo3vF43KnZlA0yHPYxL3\nORutbjNbJV76sixK1zsFi3otyJcq33a1RC++N0jG9/grx96i9oHJskb1ua9rl9C3lJldJv8AOkjY\nWwxYSoqSXvVE/puye0tRwqx16yK70YKJ2IU5rOfUPYksl7GOe6/jdhIqXL6Ev/JAODrWezemyRIk\nk7WOjRF+s5EuXVehcpKlZ1mY5ER8Mm9uuzZkXsVCi6ypZFifJIzfIGMbHdfcuC5FVV13IXaOmWBZ\n3vVFfPJvjrsyZERE9iIBaIJAAAACCQBkW8lRi6Zd8j3q9uTAW/wk47zWM+3+CJtbf1IaAEgAAAAM\nyl84bQ/Bh/eaZmUvnDaH4MP7zTAAAAAAAAAAAAAAAAAHm7U4cbrYekPN2pw43Ww9PTee/SVtVlbF\nRMa6XCuctyYKHXVzNnsiWVl+C6NVS8p7pPEQ85ew7W+bv/UZxwkxxy/0dO5vxU+tDaMXc34qfnIb\nROo9Wk4CSCTioAAAAAAAAAAAAAAAAAAAAAhVRM63BFRUyLeVLRp1qomRo1FTfGq6/ibfl/I6nVUd\nHKsEcNzEVVVE1K5bk/8As4GgQ5jXKiuai4K3pemZSjjByYSrG3vWI9bnX3ovEmTKuRSHWg9rmN3p\nqq5L7kfeudE9/wCQGgCtSPfUUSPeqK5963cSIuZDi2lejkXBiyL/AFbQFTXx073tVrnb0zfJFb9F\nvL+S+wtoqKl6GXX0U0ktWsLUd3VAkV9/gql+VfR335F5iSYD47sBGojWPvvvyZ7gOt1fC10zb1vh\nwUdyXrmO2CZJsPvVa5jsFyem6/3mZ3FUrLVsVjVhkaxuZEw0TwkzrnRVLlnUvcjJGNRzYlffGxy3\nq1Ls3tAugAAAAAAAzre4Im1t/UhoGfb/AARNrb+pDQAkAAAABmUvnDaH4MP7zTMyl84bQ/Bh/eaY\nAAAAAAAAAAAAAAAAA83anDjdbD0h5u1OHG62Hp6bzX6SrW6TxEPOXsO1vm7/ANR1bpPEQ85ewxu6\n6jed531293XYPFcddPTuenjt7VLe7Y3N+Kn1obRi7m/FT85DWlljhjV8r2sYmdXLcefqPUqzh2Ao\nJWVFQv8AB064H82bvU9SZ1Hcta9Pna9W+iKNET87zhuni+F8FHuGoblZaM9/9TWqnYRhWjT5XJHV\nN/p7x355F/Ibm/8Ai+CtTV0NQ5WIqslTwo3pc5PUWS7rvukABQAAAAAAAAAAAABAuTkJAEXJyHBI\nWJKsiJ3yojfV/wDKdgAjMSABAJAEAkAAAAAAAAAZ1v8ABE2tv6kNAz7f4Im1t/UhoASU6u0aekk3\nuRXq5G4bkYxXYLeVbsyFsx6qRKO0qySaN7mzwNSPBYrsJUwu915U9oGu1yPYjmqioqXoqcZXp66n\nqamenhfhSQXYeTIl/p9SlSlZUQ2NFQxKiVsdKxL3ouCi3XZ7vQVrIhqqe2KlklPFHEkEaYTZFdeq\nK7jVqXrlW8C3S+cNofgw/vNIx2JULuhru53xN+Zhvw2qv1+RUNdt+CmFdfx3AcgAAAAAAAAAAAAA\nAADzdqcON1sPSHm7U4cbrYenpvNfpKtbpPEQ85ewwDf3R+Ih5y9h5ySTBua1L3rmQ9XT5THS3rGV\n2bNj1aU8UrGMWWd6pgRt4/SvInpNWGhc+RJ61yTSp4LfoM1J7yhuXiayOdy5XqqXuXOpsVE7KaFZ\nZL8FORDwa0uWpd1k3nd23IDOS26NVRMJ+X+k0c5nLG48xsBIMivVUkNU1N8b3zcrXtW5zV9ClZlT\nLRPSKtXCiVbmTomTU7kX0mgcZI2yMVj2o5rkuVFTIpNks+HJFBmsc6zJEjkcr6Ry3MeuVY15F9HI\nppIuQSku6QQSVQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAIAz7e4Im1t/UhoGRbzanF8y4cO9XtyYC\n4WdOO/3GuBJBIAgEgDMpfOG0PwYf3mmZlL5w2h+DD+80wAAAAAAAAAAAAAAAAB5u1OHG62HpDzFs\nyNitffHeC3BVT0dN5r9JVrdRJgU8F2VyvVETlyHm3SR0rMOV173LxJerl5EQ2bWjkdEysqe9e9Vw\nWr/42cmvlMOjbvz3Vb0yvyR38TOL25zv00twlc73u9eh3L10EizQKqsmVUVI3pcqp6OU07b4Mk1p\n2mVZdGlZRzIjt7mjej4pEzsdcuUs1VWtXYkiyNRk0b8CVnI5F7OP1nHL1/26ezDZ4bdZ7duZDxDP\nDbrPbJmQ6dZzExcgQDxNJAAHW9jZI3Me1HNclyovGUYHuoZ20kzlWJ/iJF/SvuNI6aqnZUwuikTI\nvHxovEqEsSz4dqElGhqH4TqWp8fFx/XbxOLolJd0ggkqgAAAAAAAAAAAAAAAAAAAAAAAAAAAADOt\n7gibW39SGgZ9vcETa2/qQ0AJAAAAAZlL5w2h+DD+80zMpfOG0PwYf3mmAAAAAAAAAAAAAACCSAIv\nPKV8jau3muat8THtS/6zk2G9Xyve5tHA5UllzuT6DeNfchjV0TILWiijTBY3AREO3TzxZ2f4xl3d\nm7NVWipoUzzS4HquvX8kUyURGoiIlyIaG6B7p7VY1PFUzLv73bET8yge3pJ/Dcyb25vxU/OQ6N0L\nHUku/sT5mqwY5U5HJ4Lvd7Dv3N+Kn5yHduki3+xZ40W5VuuXkW/Ip5dT1+zU4ecZ4bdZ6W2aiWmp\nGPherXK66/1HlaKXf4onqlzszk5FTOh6bdB5DHz07D0622WpgzFuy5Xz2fHJI7Cet96+tS2UbF4L\nh9fapePBqTbOtxIAMAQSAKdfTOka2aDJURZWLy8qL6FO2kqG1UDZWZL8itXO1eNFO5TOn/gKvulM\nlPMqJMn1XcTvcpL27s3t3aJJCKCtJAAAAAAQAJBBIAAAAAAAAAAAAAAAAAAAZ1vcETa2/qQ0DPt7\ngibW39SGgBIIKlXaNPSSYEivVyNw3YDFdgt5VuzIBcB1rK1Id9RVczBwu9S+9PQdNNXRVErokSRk\njUwsCRitW7lQCtS+cNofgw/vNMzKXzhtD8GH95pASAAAAAAAAAAAAAHTUzspoHyyLc1qe30Hapn+\nX192enpnf5SfDt1EqWogVtHC+rrVwZZlRXceDyN9RkVs8dTa7JInYTFVqXmrb/B/96HlqiR0MEkj\nMr2p3qcq8R7ul05MLmze3Zo1CK+jfVLnqKp7k5qd638m/mUjXtSDuWyaGD+WiNVeVbjIPR0vpl5b\n25vxU/OQtW5wZJrTtKu5vxU/OQt25wZJrTtPJn/Y/azh42m+ZtF8P0ZbpW68zvcvrPVbofIo+f7j\ny1WisdDUNzxSIq81ci7fUep3QZaGJUzYadh31O2rjCcLNi8Fw+vtUvFGxeC4fX2qXjw6nnv2sSAD\nCgAAg4SRtljdG9qOa5LlReM5gChQSOhkdRTOVXxpexy/TZxetMxfQp2hTvkY2aDx8K4TP6uVvrO+\nlnZU07JmeC5PZ6CT4Znbs7gQSVoAAEFaprqelcjZn4Kql6ZFUsnn90flEXMXtOujhM89qluz0DVR\nzUVMy5SThD4lnNQ5nJQAAAAAAAAAAAAAAAAAAZ1vcETa2/qQ0DPt7gibW39SGgAMeqk7ktKsklje\n5s8DUjVrFdhKmFe3Jx5U9psADPs53cVm09NPhb5BTMw7mqvFd68xXoO/th8sL5Z4liudJK1UwFvy\nNbkT0+w2ABkRpMu6Gu3l7G/Mw34TFX6/pQ1m33JhZV47jOpfOGv/AAYf3mkBIAAAAAAAAAAEEnFz\nka1XOW5ES9VAq2jO+ONsUPj5lwGejld6jupYGU0DImZmpn5fSVaBrqiV9dIiphpgxIvEzl9ef2F8\nk+WZ37szdBwf/eh5WVb6iki/mVDE9i3+49Vug4P/AL0PJOcq2nBg54mq/UuRE959Hp++jZ/peXp9\n0fiIecvYYBYqa2era1Jn4SNypkK56NHC4YeGs2t7c34qfnIW7c4Mk1p2lTc34qfnIWrc4Mk1p2nh\nz/sftr2eWViSosbkva/vV9ZpzyudYcUMrr5aaXeXqudVTMvrS5fWZzPDbrQv281aavbcnzdWiOVf\n62pd+adh6db1cEnDasXguH19ql4o2LwXD6+1S8fP1fPftqJABhQAAQAABn+QV92anqnf4yfHt1mg\ndVVTsqYHRPzOTIqZ0XiVCVLHaSU7PqHyxujm8fCuBJ6V5dSlsQl3SCCSqg8/uj8oi5i9p6A8/uj8\noi5nvPR03qJlw3ofEs5qHM4Q+JZzUOZwvKgAIAAAAAAAAAAAAAAQSAMi3mz4umXfI97vb3uAt/hJ\nx3+41jPt7gibW39SGgBIAAAADMpfOG0PwYf3mmZlL5w2h+DD+80wAAAAAAAAABAEmdWr3VO2hYve\nr306pxN4k9fZeWqqobS075n3qjUyImdV4kOqz6d8MSvmyzyrhSL6eT1ZiXv2ZvfstNREaiIlyISE\nClaZm6Dg/wDvQ8hT9/XVMnJgxp6kv956e2a2nnpN7ilRz0cmQ8xZ3fU6yfzHud6r8n5XH0emlmMl\nYq0AD2st7c34qfnIW7c4Mk1p2lTc34qfnIW7c4Mk1p2ny8/7H7bnDy7PDbrQ2N17XYqjlYl74pEe\nnqTL+V5js8Nus9Buh8hj5ydh6Neb6mMScLFhuR9kwOat6KiqntUvmHuTlRLNdRr4VM9Wpf8AVVb2\n7PUbh8/PzVtIOEsrIY1kkdgtTOpwgqIqlquhej0RbluJtdtx3ArQ1tPPJvcUqOfyFgWWciQQcZJG\nxRq963Nal6qQcyFOmnqoalFWF6PRM9x3CzbkZ9ci0s7K5iLc3vZkTjby+raX2qjmoqLei5lDmo5q\ntVEVFzoUKBVppn0L8zEwoVXjZyerN7DPFZ4rQJIJNNIPP7o/KIuZ7z0B5/dH5RFzF7T0dN6iZcN6\nHxLOahzOEPiWc1DmcLyoACAAAAAAAAAAAAAAAADOt7gibW39SGgZ9vcETa2/qQ0AJAAAAAZlL5w2\nh+DD+80zMpfOG0PwYf3mmAAAAAAAAAIUFS0J3xxtih8fMuCz0cq+oVLdnUn8daF//gpl9Tn/AA7T\nQuOqlgZTU7IWZmpn5fSdxISIC5lJIXMpVeErpN6jnkX6KOU40ke9UsUf1WIhwtZMKJ8af+SRGepX\nZfyvLB9rHlzAAdEb25vxU/OQt25wZJrTtKm5vxU/OQt25wZJrTtPl5/2P23OHl2eG3Weh3QeQx89\nOw88zw260PQ7oPIY+enYenW9XBJwqUX8I6gq0T5uZFp5V9aqxfbenrPQmTTUqVm55IFXBV7VwXJ9\nFb1VF9SlyzKruuhjkcl0iXskbyOTIqe0+fq+e/bUcbX4Nm1FTc55LLz/AHFu1+DZtRU3OeSy8/3H\nWehfs91OxOFX6ndp6M85YnCr9Tu09GXqfP8AogVrS4On5ilkrWlwfPzFOGPmis7c34ufnIbRi7m/\nFT85DaOnUerUnAU7QgfLE2WC7uiFcOP08qetC4DhSzd1UtQypp2TMzOTNyeg7TPT+BtC7NT1K5P6\nZPj2mgghKHn90flEXMXtPQHn90flEXMXtPT03qQvDeh8SzmoczhD4lnNQ5nC8qAAgAAAAAAAAAAA\nAAAAAzre4Im1t/UhoGfb3BE2tv6kNAAVKq0qekk3uTDVyNw3IxiuwW8q3ZkLZj1UiUdpVkksb3Nn\ngakatarsJUwu915U9oGo+ZjYFmvVzEbhd6l96ehEznVTV8VTI6NrZGSNRHKyRisW5ePKVqGTuCyo\nYJ0eslPTMV6Naq5kuycq5DhZMqVdRJVyK5J3tREjVqokbOTLnXlA50vnDaH4MP7zTMdiVC7oa7eH\nxt+Zhvw2qv1+RUNdt+CmFdfx3ASAAAAAAgAcXvaxqucqI1EvVV4ilQNWolfXSJdhpgwov0Wcutc/\nsIrFWsqW0TPFtudOvo4m+vsL6IiJcmZDPNZ5qSSCTTQQuZSSFzKB8/ru+tGnj/re9dSJtVCwV39/\nbEn/AK2Xe1V2Fg+1h8udAAdEWaWuno2uSFUTCyrel5uWq5X2LhOzuRqqeaPR2lwE3mtPFr4yZ42f\nLUedZ4bdZ6HdB5DHz07DzzPDbrQ9Dug8hj56dhrW9XAnDnZ73R2Ej2+E1rlT2qZtiV8iWvNHMqIy\nq75qolyb4iZU9aJf6lNCi83v7H9qmZZ9K6qpKlsaok0atkidyPS9U2HmyxlmdvyrlaVoVDpp6dXJ\nvaOVt13FeX9znksvP9x599R3VI+bBVqucuE1c7Vvyp6lPQbnPJZef7jtrYyaPYnKnYnCr9Tu09Ge\ncsThV/Nd2nozh1Pn/SxJVtLg+fmKWirafB8/MU4Y+aKztzfip9aG0Yu5vxU+tDaOnUerUnCQAcVd\nFVTsqoHxPzOTIqZ0XiU6rPqHSxuimyTwrgP9PIupS2UK9jqeVtdE29WJgytT6TNqZ/aS/LN7d188\n/uk8oi5i9pvMe17Ec1UVqpeipxmDuj8oi5i9p6em9SF4b0PiWc1DmcIfEs5qHM4XloABAAAAAAAA\nAAAAAAACAM+3uCJtbf1IaBkW8lTi+ZcOLer25MBcLwk47/ca4EkEgCASAMyl84bQ/Bh/eaZmUvnD\naH4MP7zTAAAAAQAK9bU9zQK5G4T1XBY36zlzIWFVERVVciGfSJ3bU92uT5pt7YEXk43evsJUvw76\nGmWngueuFK9cKR3K5SyEJKSbAAChC5lJOL1RrFVcyIB4CDvqytk5ZcFNSJd23lgr0S4VMkn81zpP\naqqWD7Wn5Y50AB0QN2uq6eSx0iZK1X3NyIYQOWenM7L8LKlnht1nod0HkMfPTsPPM8Nus9Dug8hj\n56dhx1vVwWcJo/N7+x/apX3N55/V7yxR+b39ju1Svubzz+r3nny8mp9r8MG0P4W26jiinlX1P+J6\nTc55LLz/AHGHa0TZqqqjemRz3dppbkqhVp5qaZyd0RuTC/qS7I46a3bS2Jy42Jwq/U7tPRnm7E4V\ndqd2npDj1Pn/AEYpKtpcHz8xS0VbT4Pn5inDHzRpnbm/FTa0Noxdzfip9aG0dOo9WpOEgA4qHFyX\npcpyIUDIpqynoKp9HJM1IkVViVfo8rfVxFG3aiKomjWF6PRG3LcUrVZvlbUJmXDVUXkUrRPw2ZUu\ncmRycinv0dGYZS7uW/s9zD4lnNQ5nCHxLOahyPDeXVIAIAAAAAAAAAAAAAAAAM63uCJtbf1IaBn2\n9wRNrb+pDQAkAAAABmUvnDaH4MP7zTMyl84bQ/Bh/eaYAAACFBVrqpYI0bG3DnkXBjb6eVfQgS3Z\n01blrJ+4Y1VGJlncnEn1da9hea1GNRrURETIiIdNHTJTQ4N6ue5cJ71zudxqWCQk9wkgkqgAAFG2\nplgsiqe3wt7VE1rkQvGZbab5HS0386pZfqauEv6QPLsj3ljY/qJg+w5HKTxr+cpxPuzhzAAVAAAS\nzw26z0O6DyGPnp2HnUW5UXkLdXaU9ZEkcuBgot+RLjz6mncs8cp7LL2bFF5v/wBj+1Svubzz+r3l\nmi83v7HdqmJR101Fhbzg99nvS888wuczxny0i0OEKj8R3aXLNo3z0qz0zkZVwv8Am3LmVLsrV9Cm\ndLI6aV8j7sJ63rcWKO0Z6ONzIkZc5b1vS876unllp+GJL3c9zc6T2i/JgSNRyPYudqnqjwt8sNSl\nXSvSOoT6Spejk5HJxoXqXdNVPqI4KlscUjnIlypkdqU82vo53LdZXrCtaXB0/MUsla0uDp+Yp5cf\nNGmdub8VPrQ2jF3N+Km1ocqS06ia1HU78DARzkyJlyHfWwuWplZ7JGwDhLPFAxXzSMjamdXOREKL\nrbpVXBpmT1S/+iJXJ/lm/M8ytEGclfWvS9tlSon9cjEXtIW1Xwr/ABdn1ULeN7USRqf4rf8AkBhW\njwhPz1KEqKx2+tTnJyoW6qaKoq5pIXtexXrcrVOo+xMfFpxyym72lO5HU8bmqiorUuVDsPM2LaEk\nNUykkVrad9+B/SvIemQ+TnjcbtW5d0gAy0AAAAAAAAAAAAAAAAzre4Im1t/UhoGfb3BE2tv6kNAC\nQCnV2jT0kmBJhq5G4bsBiuwW8q3ZkAuA61lYkW+33swcK9EvvQ6qKthro3SQ4VzXK1Uc1Wqip6FA\nq0vnDaH4MP7zSM2l84bQ/Bh/eaQAAh7kY1XOVEaiXqq8QHXUTspoXSyLc1qFahge+R1ZUpdM9Lmt\n/lt5NfKcIGraE7aqRFSnYt8LF+kv119xooZ5Z57lxIBpoAAAgkgAZtSqy29SxZ2wwvlXWtzU/cXU\nqIVk3tJWK/Ng35SjQqk1q2jP9RWwJ6kwl/UB5yTxr9anE5SeNfzlOJ92cOYACoAAAAAPR0Xm8vMf\n2qecN+kqIW2GsbpWI/Ad3qrlzqYB5dCfyz+2qAA9TIS2KOZ7I5WI9quS9FT0kHZB5RFz07TOXlqv\nRYsmgX+BrpYm/wAqT5xn55U9pXtDG7KGZHJRytwVvXvmL7zaK1pcHz8xT4uPmjowLAxo6OZI+5I0\nvS9XYTthXoaSpmtlzZa17O+dfvDUb+a3qa25vxc+tCtZ3Dz+e89ln8tRlqRWPQxvSR0KTSJmfMqv\ncnrUvIiNS5ERE9BJwllZE3Cke1icqrceFpzB0d20v8+P/JDnHPFNfvUjX3Z8Fby7UeLtal//AKVR\nLTKkcuGt/wBV+tDpp6hJ0citVkjFuexc6KX7R4Qn56mbVRORyVECfOsTN9dOQ+vhPDjLGGtZMDKm\nqWKRL2uYvq9Ju0NQ9sjqSpX5+NMjv5jeJ20xNzkrJqxkjFva5iqb9bS90Ma5jsCaNcKN/Iuw+f1f\nqbw294tXklWiqu6GOa9uBNGt0jPqrsLR52pdwABQAAAAAAAAAAAABnW9wRNrb+pDQM+3uCJtbf1I\naAAx6mTuS0qySWN7mzwNSNWtV2EqYXe5Nae02SAKFl/wln01LNhb7BTMw1uW7Il2f1HTYtRHJLWN\nbhorp3PTCY5t6ZMuVDVAGIsqx7oa3+Jhgvhh8Yl9/h+lDSbXUuCl9VCq8a4aFGOngn3Q12/QxyXQ\nw3YbUW7wy66goWpetJToif8AraBPd9H9qh/zQznVtNaUytWpiSjjXLe9PnVT3J+Y7hpbTkujpYWU\nbVyvSNEWXV6PTxmilnUSJd3HB0abCcs+YSvo0S7uqD/NCcYUf2qDpEGL6L7HT9E3YMX0X2On6Juw\nrR3fR/aoekQju+j+1Q9IhOL6L7HT9E3YRi6i+yQdG3YBxo7RpK58jaaZJFj8K4toZFhUcFPLXOiY\njV7oczJyIiLd+Zrkm/u3n4fF/DhJBJBWHm6dUTdA5V4nu95p2El9m7/nWoe+a/lwlVU/K48/aEro\nqypwFukkc6Ji8jnd6naesp4W09LHCxLmxsRqJ6ES49HUc4/USPHS+NfzlOJyl8a/WpxPqzhzAAUA\nAAAAAAAAAAOcHlEXPTtOB2QeUR89O0zlxVe0Ktp8HT8xS0VbT4On5inxcfNHRn7m/Fz60K1n8PP5\n7/eWdzfi59aFazuHn8957cvNqfTPw9GZW6HyFvPQ1TK3Q+Qt56Hl0fUi15w29zfhT6k95iG3ub8K\nfUnvPo9T6VYnLNtHhCfnqVizaPCE/PUrHXT8kK7LDe2htxuEuDDUXonI2T49us9oeSs2lZWTvgkz\nPYuVM7V4lT0ob9lVL56d0dRkqYHb3KnKqfS1KmU+Z1U21G8eE1lM9z21NMqNqGJdlzPT6qnbSVTK\nqLDaitci3OY7O1eRSxcUqqmkbL3VSXJOiXOauaRORfTyKeXjul7d4uoSV6Sqjqo8Jl6Oatz2OyOa\nvIp3mmkgAAAAAAAAAAQSAMi3mVGL5l31m93t73Ay5047zWM+3uCJtbf1IaAEgAAAAMmKVkNuWjJI\n5GMSCFVVf7zmjZbTVHSI6OjTMxcjpdfIno4zoZSxT7pax8qK7e4oVRqrkvvfluNlCbbptuhrUa1E\naiIiZkQ5AFUAAAgkgDPsjwq7/lO7GmiZ1keFXf8AKd2NNEAZtrV8lDvW9ta7DvvwjSMLdJ/+f+73\nHXQxmWpJUvDz9RUufatM/Bart8WdWrmyZvzVD11lVklbTvkka1Fa67vdR4uD52snm4mfNN9WVfz7\nD1e53yKXn+5D16+OPg8TMYEnjX85TicpfGP1qcT2zhAAFQAAAAAAAAAAAljlY9rkztW8gE5Grj+p\n/lxfntNOWZ1RYb5XIiOdGqqiHlz0bPNz/qU8WvpYYeG4z3aldW5vxc+tDM7ofSWlLKxEVyPdkXWa\ne5vxU+tDGqvKpue7tNYSXVzlPZ6ay6t9bTLJI1qKjru9K+6HyJvPJ3PeQO56kbofIW8880kmvtPl\nr2ecNvc34c+pPeYht7m/Cn1J7z29T6VYnLNtHy+fnqVizaPCE/PUrHXT8kK0rB4RTmqcrYrpbNty\nOeJrd7c1sc1/Girkd6uxTjYPCKc1SbdY2W097el7XNaipyoeTUwmetZfhZw9KmVCLjPsqV8SyWfO\n5Vkp7sBy/TjXwV9y6jSPntsK3JJKOphnpsFsyouEq5npyKaFn17auNEcixzI1FdG7PrTlQzt0nhw\nalLyUbKmip3XqyVjEwJG+E3Idc8JNPHKe7HF7L94KDKx9M9Ia9EbetzZk8B+vkUvopxl3al3SACq\nAAAAAAAAzre4Im1t/UhoGfb3BE2tv6kNACQAAAAGZS+cNf8Agw/vNIzaXzhtD8GH95pgAAAAAAgk\ngDPsjwq7/lO7GmiZ1keFXf8AKd2NNEAee3WSpBTxSL9FHetch6E8punlbVWhT0zcrae98nJhL4Ke\n/wBh20N/+k2SsujiWGmYxfCuvdrXKp6fc75FLz/ch549Dud8jl5/uQ93UzbS2ZnLAl8a/WpxOUvj\nX85TiemcIAAqAAAAAAAAAAAAAAejZ5uf9SnnD0bPNz/qU8nVcY/bUV9z8sccc2G9rb1TOtxk1Sot\nVMqLeivXtOoHbHT8Ody+U3egsGaKOiVHyNauGuRVFvTRyUbUZI1y4eZFPPg5fjz/AKePddw2Nz0j\nI3T4b2tvRLr1uMcHbUw8ePhSLFeqOrp1aqKivW5UK4BrGbSQaFiPay0EV7kamCuVVOdrPY+1Wua5\nHN73KimYDndL+fjN3o7UczvK2lkY6pp71RuEnzjV8Jv/ANxoaFLUxVdNHUQOwo5EvaqHjDXs5Vsm\nCKZb+4Z8snJC/wCtzV4+RTwa3T/8pLu3Lu7N0nhQalNej8jh5jewx90aoqwKmVLlNij8jh5idhM/\nRx/ZOXY9jZGKx7Uc1UuVF4yj3NUUWWidhxfyHrm5q8Wo0Aeaws3Vaavinfva3xTJnikS53x9RaOm\nopYalmDNGj04r86alK3c9ZS+TTJNGn/jmXL6nbR3id40AUEtNjFwauKSmXlene/5JkLjJGSNRzHI\n5F40W8brLK5ggFVIIAGfb3BE2tv6kNAz7e4Im1t/UhoACpV2jT0km9yYauRuG7AYrsFvKt2ZC2Y9\nVIlJaVZJLG9zZ4GpGrWq7CVMLvcnHlT2garpY0hWZXpvaNwsK/JdynRSWhDVvVke+NcjUciSMVqq\n3lS8oLBJ8ncW3OWpZRtRUuW5Vuuuv9RyppErLWp54WvSOGnc1+ExW985W3Jl5MFQO2l84bQ/Bh/e\naRmzUFXjGWqpapkW+saxzXxYXg3+lOUnue1ft8HV/wDYDRBndz2r9vg6v/sT3Pav2+Dq/wDsBoHG\nR7Y24TluS9E9q3FDue1ft8HV/wDY4y0dqStRrq+C69FyU/It/wBYDTBndz2r9vg6v/sO57V+3wdX\n/wBgJsjwq7/lO7GmgY8Nm2nBvu92jEm+vWR38Pxr6/QdFRYtrVLvnLdkwONjIkai+xb/AMwL9baD\nkkWloWpNVqn9sfpcvFqzqZVq0bKKmpomuV7lVznyOzvct16qXqazK2kj3unqqWJt96o2mzryr3xw\nrLJr63A320Iu9vuug+J10c5hnMqlYJ6Hc75FLz/chS+TVVpBnQ/Et0ll2hRxujir4rnLet8HxPRr\n6+GeHhiSMKXxr+cpxNV25uqc5VW0GZVv8T8SPk1VaQZ0PxO06rTTwspzkbdfxrchJpO3L1LsG+0G\n96t6fM/E5fJqq0gzofiPy8DwssGp8mqrSDOh+I+TVVpBnQ/Ev5emeGssGp8mqrSDOh+I+TVVpBnQ\n/Efl6Z4aywanyaqtIM6H4j5NVWkGdD8R+XpnhrLBqfJqq0gzofiPk1VaQZ0PxH5emeGsshXIjkaq\n5VzGr8mqrSDOh+JxXcvUq9rsYNvbfd8z8Sfl6Z4azT0bPNz/AKlKHyaqtIM6H4lxLMtFKPuXu+He\n8HB8Rl/UcdbXxz22WR58Gp8mqrSDOh+I+TVVpBnQ/E7/AJemnhrLBqfJqq0gzofiPk1VaQZ0PxH5\nemeGssGp8mqrSDOh+I+TVVpBnQ/Efl6Z4aywanyaqtIM6H4j5NVWkGdD8R+XpnhrKwkw8G/LdeSa\nXyXqcPDxg2+67xPxOXyaqtIM6H4k/L0zwss9XZjWvsqJj0RzXNVFReNDI+TVVpBnQ/EvwUNp08LY\nmV8OC3Il9P8A7Hn6jWx1MZMVk2Y1uwS2Y6JrFfNSZcFl174vQnKnoznorLqoKqhhdBK19zERblyo\nt2ZU4ihWWRX1qtWW0Iu9zYMF3vOuXc/USoi91wMkRqNSWOBWvu1o44ZZS4TH4VvgwqeyLZgyY+fI\nnJJC1fzzlttNayZ7RhXXT/E5q0gZ3c9q/b4Or/7E9z2r9vg6v/sBfVEVLlRFT0lF1n0j5nJEjoJU\nuVVhcrM+rIce57V+3wdX/wBjilHaiSukSvgwnIiL/D8n93pJtulkrt7nr4Uugq2yJyTsvX2pcO6a\n6PxtEj/TDIi/ktxx7ntX7fB1f/Yjue1ft8HV/wDYbJs540iat0sVREv9US3e1L0Oxtp0LsndUSLy\nK67tOjue1ft8HV/9jg+htB/h1dK7XTX+8dzum25Y5LIlwHtdlbmW/wCkhpmDNYEs6XPlpEyot7KV\nGuyLfnvN4qzf3ASAqASAIBIAgEgDhJIyKN0ki4LWpeq8hVjtSikkaxk6K5y3Ily5S4LgMOe0JJbX\npFiqEZTb86JWoqfOKjVvVfQipd7TcKM1kUc1RDMsMbXxPw0uYnfLcqZfaXgJAAAgkAQCQBB1zzx0\n8SyTOwWJnW47SAKcNp0dTIkUM+E9cyYKp2oVmb5Da0MEVTNP3quqEeqKjUuyL6Fv4tZqql6chQob\nNfRSK5KyWRrnK57XNb3yryrdeBfBIAgEgCASAIBIAgxK100r7SlbUSRrRx/NI1bm3o3CvXlNso1V\nlxVUsj1klYkrUZKxq5JETiX4AWqaRZqaKRUuV7Ecqa0OwNRGoiJkRMiEgQCQBAJAEAkAQFJIAyqr\nfYK6nbBVTSTSS3uicqK1I7++yXZETiXluNUz4rMfDWyVLK2W+V+E5qtat6fVvuvuNEAQSAIBIAAA\nAAAIVbkvUxIK+Sotqnc2dEpZWPRkSKnfXXXOXXl9Rtql6XLmKWKaNKyGpZDGx8SKiI1iJnuy/kBe\nIJAEAkAQCQBn49sfStD1hm0Y9sfStD1hm0+HgD7hj2x9K0PWGbRj2x9K0PWGbT4eAPuGPbH0rQ9Y\nZtGPbH0rQ9YZtPh4A+4Y9sfStD1hm0Y9sfStD1hm0+HgD7hj2x9K0PWGbRj2x9K0PWGbT4eAPuGP\nbH0rQ9YZtGPbH0rQ9YZtPh4A+4Y9sfStD1hm0Y9sfStD1hm0+HgD7hj2x9K0PWGbRj2x9K0PWGbT\n4eAPuGPbH0rQ9YZtGPbH0rQ9YZtPh4A+4Y9sfStD1hm0Y9sfStD1hm0+HgD7hj2x9K0PWGbRj2x9\nK0PWGbT4eAPuGPbH0rQ9YZtGPbH0rQ9YZtPh4A+4Y9sfStD1hm0Y9sfStD1hm0+HgD7hj2x9K0PW\nGbRj2x9K0PWGbT4eAPuGPbH0rQ9YZtGPbH0rQ9YZtPh4A+4Y9sfStD1hm0Y9sfStD1hm0+HgD7hj\n2x9K0PWGbRj2x9K0PWGbT4eAPuGPbH0rQ9YZtGPbH0rQ9YZtPh4A+4Y9sfStD1hm0Y9sfStD1hm0\n+HgD7hj2x9K0PWGbRj2x9K0PWGbT4eAPuGPbH0rQ9YZtGPbH0rQ9YZtPh4A+4Y9sfStD1hm0Y9sf\nStD1hm0+HgD7hj2x9K0PWGbRj2x9K0PWGbT4eAPuGPbH0rQ9YZtGPbH0rQ9YZtPh4A+4Y9sfStD1\nhm0Y9sfStD1hm0+HgD7hj2x9K0PWGbRj2x9K0PWGbT4eAPuGPbH0rQ9YZtGPbH0rQ9YZtPh4A+4Y\n9sfStD1hm0Y9sfStD1hm0+HgD7hj2x9K0PWGbRj2x9K0PWGbT4eAPuGPbH0rQ9YZtGPbH0rQ9YZt\nPh4AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAA/9k=\n",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/M_kyGLaqqgQ\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7fcdf853ea90>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "YouTubeVideo(\"M_kyGLaqqgQ\")"
   ]
  },
  {
   "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
}