{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Metropolis Hastings\n", "\n", "\n", "This is a simple notebook to play with the Metropolis Hasting (MH) algorithm. We will be simulating the following distribution\n", "$$ p(x)=6x(1-x)$$ with $0\\le x\\le1$. This example is based on notes linked in reading.\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAECCAYAAAAYfWtSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X18zfX/x/HHOds529l2tmFzleurNyKlIhcpoiiiEKLS\nBSWRQkVERYpfQnKZUunqmy6EJFJCKqUS6W1Gri/Grs7O2Xa2nfP7YxsjbDu7+Ozidb/durV9zufi\nubftvM7n835/3h+T1+tFCCFE+WQ2OoAQQgjjSBEQQohyTIqAEEKUY1IEhBCiHJMiIIQQ5ZgUASGE\nKMf887KSUqo18LLWuqNSKhJYDIQDfsC9Wuv9SqkhwFAgDZiqtV6tlAoElgGVgUTgPq316aL4QYQQ\nQuRfrmcCSqmxZL7pB2Qtmg4s01rfCEwEGiulqgAjgDZAV2CaUsoCDAN2aK07AO9lrS+EEKKEyMvl\noL3AHTm+bwfUUEqtA+4GvgdaAZu11ula60QgCmgBtAe+ztpuDdC5kHILIYQoBLkWAa3150B6jkV1\ngFitdRfgEPAMEAok5FgnCQgD7DmWO7LWE0IIUUL40jF8GliZ9fVK4Boy3+hzvsHbgTgy+wHsOZbF\n+xZTCCFEUchTx/B5NgG3Au8DHYCdwDZgqlLKCtiAxlnLf8xa99es/2/KywG8Xq/XZDL5EE0IIcq1\nfL9x+lIExgBvKqWGkXkGcLfWOkEpNQfYnBVivNbarZSaD7yjlNoEpJLZh5Ark8lETIzDh2hlT2Sk\nXdoii7TFWdIWZ0lbnBUZac99pfOYSugsol75R80kv+BnSVucJW1xlrTFWZGR9nyfCcjNYkIIUY5J\nERBCiHJMioAQQpRjUgSEEKIckyIghBDlmBQBIYQox6QICCGKzNy5sxgx4mEGDuxD797dGTnyEZ57\nblyRHGvSpPH88cf2C74WFbWHpUvfBOCHH77n9OlT/1nn9ddncvLkCd56axErVnyWp2O63W5WrfoC\ngDVrVrFlS57uhy1RfLlZTAgh8uSxx0YBmW+QBw8e4OGHhxuSo2HDRjRs2AiATz75kDp1xlOpUsQ5\n64wY8WS+93v69ClWrlxB9+696Nate6FkLW5SBIQoJ4InTyBg5ReFus/UHr1wTp6S7+1+//035s9/\nHavVSo8evXjzzQV88MGnWCwWFiyYS+3adejWrTsLF77Bjh1/4PFkcNddd9Ox47kTEX/66f9Yu3YV\nYWEViI+PA+DQoYO89NLz+Pv74/V6mTRpCocPH+KLLz6la9dbiYraw5Qpk5g48QWefXYs4eEVuO66\ntmzduoWxY8cDsHHjd2zYsI7U1FRGjRpD48ZN6dnzFlasWAtknnXccUcf1q5dw4ED+1m69E08Hg+V\nKkXQs+edzJ07ix07/sBkMtGlyy306dOfl156HovFwrFjx4iNPc2zz06iYUNVwH+BgpPLQUIIQ6Sl\nuZk7dxG33HIrF5ry5qeffuTo0SO88cZiZs9ewLvvvoXTmXTm9bi4WJYv/4hPPvmEadNeJS0tc7Lj\nbdt+pmnTZsyaNY8HHhhKUlLmNiaTiTZt2tOwYSMmTnwBi8VCXFwcr732BnfffS855yurXv0yZs+e\nz9NPT2D69Jeylv434333PUCdOvUYPPihM8t+/HEzx48fZdGipbzxxmLWrVvLvn17AahatTozZ75O\n7953sWLF5wVtwkIhZwJClBPOyVN8+tReVGrVqp3ju7PT12RPZbNv3160/oeRIx/B6/WSkZHBsWPH\naNCgIQBHjhymXr36+Pv74+/vT5MmTQHo3r0n77//Dk8+OQK7PYShQx/9z7Gzj1GtWnX8/Pz+8/qV\nV14FQN269YiLi/1PxnO/Pte//+7niisyt/f396dp02bs378fgEaNMj/5V65chb/++vOi+yhOciYg\nhDCEyXT27ScgIIDTp0/h9XqJitoDQO3adbn66muYM2cBc+YsoFOnLlx2WY0z29SoUYv9+/fhdrvJ\nyMhgzx4NwKZNG2nR4ipmz57HjTfexPvvv3vOcc1m85kikPPTf8551Hbv3gVAdPReqlSpCkBGRgYp\nKSmkpaWxf/++M9t7PJ5z9l+3bl127PgdgPT0dHbu/JNatWr953glhZwJCCEMN2DAPYwZM5Jq1aoT\nGpr5aJJ27a5n+/ZfGT58CMnJyXTocCM2m+3MNuHh4QwaNJh+/foREhJ25rXGjZswdepkLBYLHo+H\nkSOfPHNJCKBZsyuYMuU5xo4df86bcs6vjx49yuOPDyMtLY2nnsrsJ+jbdwAPPzyY6tUvo2rV6gBU\nqFCR9PQ0FiyYS0BA5hN427Rpz/btv/HIIw+Qnp5Op05dSsS1/4uRWURLOJkh8Sxpi7OkLc6StjhL\nZhEVQgiRL3I5SAivF/PxY/jti8Z87CjmEycwnziOOS4Wk9OJKcmBKTn53E0CAvEGB+MNCcETFoan\nStUz/2XUqYunZi3wlz8vUfLJb6koFbxeLwkJCSQmFuy03+R0EvD3TgJ27SRg519Y9T9YDuzH7HJd\n+vgmE2RfM/Z6MeVyGdVrsZBRpy4ZqgnpV7Qg7YorSW9xFd5KlQqUX4jCJkVAlAoORyJrtx7C483f\nr6x/SjJVdv5G1Z3bqbprOxX3acyejDOvpwUEEl+tJonVa5FYrSbOiCokV4gguUIEKaHhpNmCSLMF\n4bFYSU520aV1A0JDw8DtzjxDcDoxx8dlnjmcOIH52FH89u/Db99e/PbuxT9qDwGrVpw5XnojRdp1\n7Uhr2w53h454IyIuFFuIYiNFQJQaQUHBeLDmup7t5DGqb91A9Z+/J/LPX/BLcwPg8fMntnFzTje9\niriGlxPXoClJl9UG88W7xvyy/gPOngkAWK14K1bCW7FS5qWf5i3+u7HXi/noEfx3/In/n79j+W0b\nlm2/4P/uW9jefQuvyUR6y2tw39yV1FtuJaNJ03OPIUQxkCIgyoSAuFPU+GEttb7/iohdZycRi6/X\nmGOtOnDyqus43eRKMgJtl9hLITOZ8FxWA/dlNXB3uy1zWVoa/jv+wLJlE9Zv12H55Scsv20jeNqL\npKvGpPbqTeodvcmo16DQ43i9XhyOxELdp90eWiLHvou8y9MQUaVUa+BlrXXHHMvuBh7TWrfN+n4I\nMBRIA6ZqrVcrpQKBZUBlIBG4T2t9Og+5ZIhoFhn+likxMYE/9sWecyZgykin6rZN1F3zKdV+/h6z\nJwOvycTJFq053OEWjrW6geTK1Qotg8vpoH3zapmXgwqJKT4O64b1BKz6Euu6rzGlpgKQ1uo6kgfd\nR+rtd0BQ0H+28+X3IjExgXU/78UWFFwo2ZNdzrOXxy7i999/Y+TIR5g8+SVuuqnLmeX33dcfpZow\nfvwk3G43ixfP5++/d2IymQgKCmLMmHFUrlyFxx4bSnx8PMuW/e/Mths3bmDChKf55JOVVK1alQMH\nNLNmzSE9PZ2UlBRuvbUHd9zR56KZfv55K99++w3jx09iwoSnmDJl+gXXO3HiOHv3RtGu3fX/2f7k\nyRNce21rJk0az8KFb+epvX744Xsuv7wZJpOJpUvf5Mknn87TdvnhyxDRXM8ElFJjgXuApBzLrgIe\nyPF9FWAE0BIIAjYrpb4BhgE7tNYvKKX6AROBUfkNKURO1oQ46q3+mPqrPibo1HEA4ho05d8uPTnc\noSsplSobnDDvvOEVSL2zL6l39sXkSMT69VcE/u9DLD98T+gvP+F59mlS+g0gecgwPHXrFfh4tqBg\ngoLthZA872rXrsO3335zpgjs27eXlJSUM6/PmfMqtWvXZfjwx4HMN8tJk8Yxf/5bZ84y9u6NOjNd\nxLffrjtzs9bRo0eYOnUq06fPITw8nNTUVB5/fBiXXVaDVq2uyzXbxQoAwPbtv3LgwL//KQKtW7cB\n4PjxY/k6C8qevbRWrdpFUgB8lZfLQXuBO4D3AJRSlYApwOPA4qx1WgGbtdbpQKJSKgpoAbQHXsla\nZw2ZRUAIn9gP7qf+Zx9QZ/0K/NyppAUFs7d7f/Z360N8w8uL/PhFcTnlP27pBrd0w//wIezLPyb0\nk48IenMhtiWLcHa+mYQHhmDp1CX3/ZQg9es35NChg7hcToKCglm7dg0339yNEyeOk56ezqZNGxkz\n5uwzBjp0uJErr2x55vvOnW9m3bqvadCgIUlJSbjdqVTKGmW1du1X9OrVi/DwcCBz+omZM1/HZjv3\n7OnAgX+ZNu0FbDYbgYGB2O2ZdyVnzwz62Wef8PXXq/HzM9O48eWMGPEEy5YtJTU1lebNW/DRR8uo\nUKEiDkciN910M4cPH6JXr97ExcUybtxoYmNjadu2Pffd9yAvvfQ8nTvfQqtW15056+jY8aZzZi+d\nMmUSCxe+zbZtP7F48QICAgIICwtj3Ljn2LNH8/7772CxWDh69Cg33dSFe+99gKKSaxHQWn+ulKoN\noJQyA28CTwKpOVYLBRJyfJ8EhAH2HMsdWesJkS9+u/+m8vSp1PtqFSavl6SqNYi64x7+vaU36YV0\naSMvkl1ONm6PJbxicQzz9IebBmK6oR+1t27g8i/fJ2LdWkLWrSX52tYwbSo0v7bUdCTfeGMnNm78\njm7durN79y4GDRrMiRPHSUiI/8+8/sCZqSNMJhPt2nVgypRJDBs2gu+/X0/Hjp35/PPlAJw6FcM1\n11x5zrZBF/ideOON2QwZMoyrr76W999/hwMH/s16JbP91qxZxejRz9C4cRO++OJTAAYNGszBgwdo\n1+56PvpoGTff3JX27W9gzZpVZ84AUlKSmTjxRQIDAxk+fAjt2nW44M+fPXvpU089i8ViObP99OnT\nWLBgCZUqRbB8+UcsXbqEtm3bc+LEcd5992NSU1Pp1aursUXgPC2BBsB8wAY0UUrNBL7j3Dd4OxBH\nZj+APcey+AKlFeWKn/6H4FemnhliGdegMbsHDONI25vgAjM/FodAW1CxX0452bU3J2+5k0q7ttPo\n/fnU2LYFOncm/JpWOMdNJO36G4o1T35lzqnflRkzplGtWnVatLjqzGRtYWHhOBz/7dv45puv6dQp\n89kBAQEBNGqk2LlzB5s2beT556fx2WefAFC1ajWOHTtG8+Znt927Nwqv13POfD2HDh04M8to8+Yt\nchSBzBzjxj3HRx8t49ixozRrdsV/JoUDqFmz9n+W1a/fiKCsPpsmTZpy6NCBc14/v8815/fx8fEE\nBwefKYItWlzFokXzaNu2PfXqNcBkMhEYGEhAQOB/jluY8lMETFrrX4HmAFlnBx9qrZ/M6hOYopSy\nklkcGgM7gR+BW4Ffs/6f52evRUYW7x9aSVbu2uLwYZg8Gd5+GzweaNUK5+jRbAhXBNvDMKo1kp1W\nzGYL9pCi/aO8GPd17fi5WXPCLbGEzJqFZcUKwnv3gK5d4eWXocUFhqnmYLV6CAmOJbiQ8ptxExFh\nJyzs4v8i4eFBBAZauOIKRUaGm5UrP2X06NEcPHiQwEAL1apV4MYbO/D1119wzz33ALBmzRpWrPiE\ngQP74u9vpmLFYHr37sUnn3xMZGQlataMxGLxo1KlYPr1681jjz1Gt27dqFixIk6nk1mzXmH48OHn\n/N00bqw4eDCK66+/nkOHogkMtBAZacdsNhMZaWfRotW88spLWK1WHnzwQY4ciSY01EZgoD+Rkfas\n44UQGWnHbg8kKMhKxYrBHDr0LyEh/lgsFqKi/mHw4HvYtesP3O4kIiPtHDmy/8yxAgIshIdnXo6y\nWPxo2LAmqanJmEypRERE8NVXu2jUqMGZNsvObzabivQ9ID9F4KLDiLTWJ5RSc4DNZJ5fjddau5VS\n84F3lFKbyLx8dHdeDyYjYjIZPTqoWK6DZzG5XIQvfIPwNxdiTk3F3aARp8c+g6tTZxxJDrwxaTiS\nUnLfURFxOt2YzRkE2IzL4HKmcqJ5A0K++IK4dRsJfnES1q+/xrt2Lal9++Oc+DyerKmPz5eY6CAm\nJpYkZ+oFX8+vZJeTU6cq4nZf/D6L+HgXKSlpxMQ46NChE2vXriEoqCIJCfrM8oceGs7rr79Gnz59\nAROhoaE8//zLxMQ4SE/3EBvrpGHD5vz669OMHz/5zPLTp51UrVqVsWPH8sgjj+Ln54fL5aJHj140\naXLVOX83Dz00nKlTJ7NgwSLCwytgtVqJiXHg8XiJiXFQvXot7rqrH0FBwURGVqZ69Xq43TBv3nxq\n1Kh3JkdIiAOHIwWXy531vZ1HHx1BfHwcnTvfTGhoZTp3vo1p017g008/p2bNWmd+TqUuZ/ToMYwd\nO560tAxiYhyMHj2Ohx8ehtlsxm638+yzk4mO3ktqavqZ/NkZ88KXYiGziJZwRheBwh5WeEFeL7V/\nXM81S+cQcuoEroqR/D5gKNEdb8Prl/k5JfbUCSIrVybAZly30qmTxzCb/agYYdzoo+xhqvXr18j8\nvfB6sXz3LSEvPIf/3zvxhNhxjX6a5CGPgPXcG+vK6n0CRv+NlCRFMkRUiKIcVhhy5F9azn6eKn/8\nRIbFwu4BD7O7/xAybMHkvK3LleOxgiIHk4m0Tp2Ju6Ejge8tJXjaC4Q8P4HAD94lacYs0tq2z7Gq\nqVDvcRBlg0wlLQxhSk+j8YeLuHloT6r88RNHW9/A2kUr2Xn/KDJsxTfip8zw8yNl8IPEbt1O8uAH\n8dsbRXivWwkZ/TimBBmPIS5OioAoduF7/6bzY31p/vZruO2h/DhhFltemI/zsv+OvhD5461YiaTp\nrxH/1XrSmzTF9t7bVGjfCuvXXxkdTZRQUgREsTFlpNPk/fncNKIf4fs0+7r1Ye3iVRzpcEupGe9e\nWqRffS1x637AOW4i5vg4wu7tT8jjj2Iqpk5+UXpIERDFIuTwfjo+MZBm78whpUIlfnjpTX574kXS\n7HKNushYrbieGEvc+k2kXXEltg+XUeHGtli25HmktigHpAiIIld73Rd0ebQPlf7ZwYGO3flm4QpO\nXNPO6FjlRoZqTPxX63E++RTmo0cIu7M7QS+/COnpRkcTJYAUAVFk/JKdXDv9GVrNGIfXbOancf/H\nL+NmyKd/I1ituJ6ZQPzKtXhq1iJ45gzC7uyO+egRo5MJg0kREEUi9N8oOg/vS531K4hVzVk3/zMO\ndbzN6FjlXvo1rYj7dhMpt9+B9acfqdCxLZYN64yOJQwkRUAUuhob13DTyP6EHt6P7j2YDTOX4axW\n0+hYIos3LBzH4qU4ZszC5HIRNqAPQa++kjlFhyh3pAiIQmPKSOeKRdNpM/VJvCb4ceJsdjz8NF5L\n7o+EFMXMZCLlvgeIX/UNnstqEPzKVELvGyD3FJRDUgREobA4Erh+/FDU8rdJrFGXb1//H0euv9no\nWCIX6S2uIm7dD7g7dCRg7RrCb+mIX3SU0bFEMZIiIAos5PB+bnq8P1V+38rR6zry7ev/w1GrvtGx\nRB55K1Ui4ePPcA1/HP990YR3vQnLD98bHUsUEykCokAif/+Jmx4fgP3wv/xz14NsmfQ66cEhRscS\n+eXnh3PSiyTOmY/J5SSs3x0ELl1idCpRDGQCOeGz2uu+4JqZE8Fk4pcxL3Hg5juMjlTmZc8EmpCQ\nQGJiEcyceWt3HJUrU3XYEOxPPUH63j3EjnkGzOd+XiwJs4eKwiFFQOSf10vjjxbR/O1ZuENC2TJ5\nLqeuuNboVOVC9iMuo2PSC+25AP8RUJOQaUvo/MLjVFg0n7jd0WwZ8RyerA7+ZJeTLq0byIykZYQU\nAZEvpox0rpo7hfqrP8ZZuRqbpi7CUbuB0bHKlUBbEMEhoXgouofbeOo15vs5H9HuueHU2/QNIYnx\n/DjpddJC5DHhZY30CYg8M7vdXDf1Seqv/pj4eo3ZMPsjKQBlmDu0AhtfeYvD7TpT+c9fuGHsYALi\nThkdSxQyKQIiT/ySXbR7bhg1Nq/jZItWfPfqe6RUMu4JW6J4eAIC2TphFtG39aNC9G46PjmI4Jjj\nRscShUiKgMiVNSmRDuMepOr2Hzl6XUc2TVkoI4DKEz8/to+cxO5+Q7AfOUC3cQ9hid5rdCpRSKQI\niEsyx8Zy83OPEvH3Hxzo1IMfn5uNJyDQ6FiiuJlM7HzwSXY8NJrg0ye5bEAf/P7eZXQqUQjy1DGs\nlGoNvKy17qiUuhKYA6QDqcC9WusYpdQQYCiQBkzVWq9WSgUCy4DKQCJwn9b6dFH8IKLwmWJiqH5P\nPwL27yH6tn5sH/Hcf4YKivJF3/UQTj9/2ix8hfA7byN++UoymjU3OpYogFz/opVSY4HFQEDWolnA\ncK11J+Bz4GmlVBVgBNAG6ApMU0pZgGHADq11B+A9YGLh/wiiKJhOnCD8ztsI0P+w+9a+bB85SQqA\nAGBP196cfGkGprg4wnt3x3/HH0ZHEgWQl7/qvUDOu4D6aa3/yvraH0gBWgGbtdbpWutEIApoAbQH\nvs5adw3QuVBSiyJlOnmS8Dtvw1//Q/zgB/nloTHy+EdxDsdd/XHMnocpPp6w3rdLISjFci0CWuvP\nybz0k/39CQClVFtgOPAaEAok5NgsCQgD7DmWO7LWEyWY6fRpwvv0wD9qD65hIzj97CQpAOKCUvsP\nxPH6AkyJCYT17Sl9BKWUTzeLKaX6AeOAW7XWp5VSiZz7Bm8H4sjsB7DnWJbneWojI+25r1ROFFtb\nxMXBgDvgn90wciRBs2YRkZhISHAcwSHGdgYnOzPvVrUbmCPZacVstpSIDGBcW5hxExFhJyzMDsOH\nQpAF0wMPULHv7bBxIzRpUuyZ5P3Cd/kuAkqpQWR2AN+otc5+U/8FmKKUsgI2oDGwE/gRuBX4Nev/\neX7CdUxMEcyLUgpFRtqLpS1MjkTC+tyO5Y8/SL73AZKefRFOJZGY6CDJmVqkd6fmhdPpxm634Egy\nLofT6cZsziDAZnyGiEgMawuXM5VTpxy43VkXErr3IXB6AvanniCjYyfiV6zBU6/4ZpEtrr+R0sCX\nYpivnj6llBmYDYQAnyulNiilJmVdIpoDbAbWA+O11m5gPtBMKbUJeAh4Pt8JRdFLTiZ0UD8sv28n\npf9AkqbPlEtAIl9SBj9I0ovT8DtxnPC+PeXZxaVIns4EtNYHgLZZ31a6yDpLgCXnLUsG7ipIQFHE\n0tMJHToY69YtpPboheO1uTIKSFxS9kym50scMIjU2NNUeu3/sPfpyZEPl+OpUKHIcshMpoVDJpAr\nzzwe7KOGE7B2De4OHUmctxj8/IxOJUq47JlMwyte4PPg9X25du9hmq78iJCBd/PN5Lmk24KKJIPM\nZFo4pAiUY8GTniXwfx+S1vJqEpa+DwEBuW8kBJkzmQYFX/j6867hEwlKdlFn/Zd0njGOTVMWyHOm\nSzA57y+nbPPnErTwDdIbKRI+WA4hMheQKCRmM78+OYWj13Wkyu9buXbmBPB6jU4lLkLOBEowr9db\nJE+QCl79JSGTxpNepQpH3nyHdH9/SEy44LoORyLI36/IJ6+/hZ/Gv8oNT99P7W9X4oqsxs4HnjA6\nlrgAKQIlmMORyNqth/B4C++fqcqu7XSZPAq3LZivn36VuNNmOH3souvHnjpBUHAoQSEyDlvkT0ag\njc0vzKfTqAE0+WgRrsiq7OsxwOhY4jxSBEq4oKBgPBTO9VT7wWg6vfwUJq+XrZPmkNrsanLrsnM5\nkwrl2KJ8codVYNPURXQadTct35hCckRVjrXpaHQskYP0CZQT1oQ42k8chjUpkW2jp3CyZdvcNxKi\nEDir12Lzi/PJsFi5btoYwqJ3Gx1J5CBFoBwwu920m/wYIccO8ffAYRzs3NPoSKKciVPN+eXpV/BP\ncdF+4qMEnj5pdCSRRYpAWef1cs1rE4jYtZ2DN3Rj170jjE4kyqkj7W9mx4NPEnTqOO2eexS/ZJfR\nkQRSBMq8xh8upPa3KzndpAXbxrwk00EIQ+m7HmL/LXdSMWoXrWY8Ax6P0ZHKPSkCZVj1H7+l+dLZ\nOCtXY8vkufJYSGE8k4nfRk7i5BXXUmPzOpp8MN/oROWeFIEyKvTfKFq98hTpAYH8OHkuqRUijI4k\nBABei5WtE2bhrFKdZu/OpfqW9UZHKtekCJRBlsR42k1+DEuyi22jpxLfoKnRkYQ4hzu8Ilsmv0F6\ngI1W058mdP8eoyOVW1IEypqMDK6bNoaQowfZPeBhDt94q9GJhLighPqN2Tb2JSzJrswPLYl5fuaU\nKERSBMqYZu/MoepvWzja+gZ23jfS6DhCXNLhDl35e8DDhBw7ROuXn5KOYgNIEShDqm9ZT5OPFpFU\nvRa/PD1dngsgSoVd947g2DXXU+3XTTRdNs/oOOWOvEuUESGH99NqxrjMjuDn5pAWEpr7RkKUBH5+\n/PzMdJxVLuPyZW9Q9efvjU5UrkgRKAP8kl20fX4kFlcSv416gYR6yuhIQuRLWmg4Pz43mwxrAK1f\neZrgoweNjlRuSBEo7bxerp4zmbADe4nqOYiDN/UwOpEQPolveDm/jZyENSmRNi+OwuxONTpSuSBF\noJSr8/WnmXcEqyv4c+hYo+MIUSAHbr6DfV17UyF6Ny0WvGx0nHIhT1NJK6VaAy9rrTsqpeoDSwEP\nsFNrPTxrnSHAUCANmKq1Xq2UCgSWAZWBROA+rfXpwv8xyqewfZqWb0zBbQ/jpwkz5RF+okz4ffgE\nKuq/aLDqI2KuuFaGORexXM8ElFJjgcVA9gNoZwLjtdY3AGalVE+lVBVgBNAG6ApMU0pZgGHADq11\nB+A9YGIR/Azlkr/LSZspo/Bzp/LLmGm4qlxmdCQhCoUnIJCtE2aRZgvimtcmEnJ4v9GRyrS8XA7a\nC9yR4/urtdabsr5eA3QBWgGbtdbpWutEIApoAbQHvs6xbudCSV3eeb20nPM89sP/ovvcLw/pEGVO\nUs26/DbqBSzJLtpMeVL6B4pQrkVAa/05kJ5jUc5pKB1AKGAHcj6kNgkIO2959rqigGqvW0HtDZkz\ng/4lz20VZdShjrcRfetdhO/7hysW/5/RccosXx4vmfOWPjsQT+b1/tDzlsdlLbeft26eREbKM22t\nVg/si8Uecnb2z6CD+2j5xoukBdvZ+eLrhIQXbTslO62YzZZzMhgh2ZnZ32FkjpLQFtkZwLi2KM52\niB4zmcon0krGAAAgAElEQVS7f6fhimUktruBmPaZFxPMuImIsBMWlvn7L+8XvvOlCGxXSnXQWv8A\ndAM2ANuAqUopK2ADGgM7gR+BW4Ffs/6/6cK7/K+YGIcP0cqWxMTMNnAkpQCZTwhrPWEE/skuto5/\nlZOhkZD1WlFxOt2YzRkE2Ir2OHnJYbdbzrSFURmMbovsDBGRGNYWxdsOJn585v/o/Fhfmk0ZyzcL\nviAlogouZyqnTjlwu81ERtrl/SKLL8XQlyGiY4AXlFJbAAuwXGt9ApgDbAbWk9lx7AbmA82UUpuA\nh4DnfTieyNJs6Swq7P2b/bfcKSMmRLmRWLcRfz78NAGJ8bR+5WnIyDA6UpmSpzMBrfUBoG3W11HA\njRdYZwmw5LxlycBdBU4pqPLrFtTyt0msUZffH33W6DhCFKvoHgOo8tuPXLb1W9Qnb/F7j/5GRyoz\n5GaxUsCaGMe1/zcOj7+Fn8fNIMMWZHQkIYqXycSvo18kuWIkzd59nYr7tNGJygwpAiWd18vVsyZj\ni41h570jiG94udGJhDCEO7QC28a8hDk9jetfew5TSrLRkcoEKQIlXK11K6mx+Rtiml2N7vuA0XGE\nMNSJa9oT1XMQ4Yf3U3H6NKPjlAlSBEow/0MHuXLey6QFBfPLU6+An5/RkYQw3I6HRhNfoy7h776N\nZcM6o+OUelIESiqPh8pPPYHF5eT34RNwVZVpIYSAzGklNj3xAl6LBfuoxyAuzuhIpZoUgRLKtmge\ntm2/cKRdJw507ml0HCFKlNh6itgRo/A7fgxGymNUC0KKQAnkF7WH4JdeIKNCRX4fOQFMptw3EqKc\niR/6KGlXtYRly7CuXml0nFJLikBJk56OfcTDmFJSiJnyMqkVKhmdSIiSyd8fx+sLISAA+9jHMZ06\nZXSiUkmKQAlje2M2lu2/kXJnX5y3dDM6jhAlWkYjBVOnYj51CvtTT4DXa3SkUkeKQAni989ugqe/\nREaVqiRNm2F0HCFKh1GjSGvdhoBVKwj48nOj05Q6UgRKivR07I8Pw5SWRtKMWXgrVDQ6kRClg58f\njtlv4A0MJGTcGLkslE9SBEoI28J5WH7fTsqdfXF3lcnhhMiPjHoNcD4zEfOpU4RMeMroOKWKFIES\nwC86iuBXpuCJiCBp6nSj4whRKiU//ChpLa8m8LPlWNesNjpOqSFFwGgeD/ZRj2FKScHx8qt4K8lo\nICF84ueHY9Y8vFYrIU89gSlebiLLCykCBgtcugTLz1tJve123D16GR1HiFIto3ETXKOfxu/EcYIn\nTzA6TqkgRcBA5qNHCJ4yGU9oGEkv/5/cFCZEIXA9Nor0ps2wffAeli15fphhuSVFwCheLyHPjMGc\n5MA56UU8VaoanUiIssFiwTFzDl6TiZDRIyFZppy+FCkCBrGuXknA16txt2lHysB7jY4jRJmS3vIa\nkoc8gv++aIJmyT03lyJFwACmhHhCxo3Ba7WS9OocMMs/gxCFzfnMRDJq1CTo9Vn4/b3L6Dgllrz7\nGCB4yvP4nTiO64mxZDRoaHQcIcqmkBCSXnkVU3o69tEjweMxOlGJlKcHzZ9PKeUPvAPUAdKBIUAG\nsBTwADu11sOz1h0CDAXSgKla63I9gNf/t20EvvsW6aoxrhFPGB1HiFLJ6/XicCQCYLV6SEx0XHjF\n1m3wu7U7IV+tgsXzSRwwqFBz2O2hmEr5gA6figBwK+CntW6nlOoMvARYgPFa601KqflKqZ7AT8AI\noCUQBGxWSn2jtU4rjPClTno6IWOfwOT1kjT9NbBajU4kRKmU7HKycXss4RUrERIcS5Iz9aLr2noP\no9f33xH68ktsqHEVKeGFMyVLsstJl9YNCA0NK5T9GcXXIrAH8FdKmYAwMj/lt9ZaZ4/HWgPcTOZZ\nwWatdTqQqJSKAq4AfitY7NLJtmQhlp07SOk/kLQ27YyOI0SpFmgLIijYTnBIIB5SLr5isJ2d9z9B\nyzem0Pr9eWx76pXiC1kK+NonkATUBf4BFgJzgJznRA4gFLADCedtV7rLpo/Mx44S9PJUPBUqkPTc\ni0bHEaJcie7en9iGl1Nn/ZdE/vGz0XFKFF/PBJ4AvtZaP6uUugz4Hsh5bcMOxAOJZBaD85fnKjLS\n7mO0EurRCeBMgtlvEtGkbp42sVo9sC8We0hgEYe7uGSnFbPZYmiG7BxAuW+L7AxgXFuUpHbIzpCX\nLP+Mm0abh3pxzdwX2PLuV3itAQXKYMZNRISdsLDS/V7laxGIJfMSEGS+qfsDvyulbtBabwS6ARuA\nbcBUpZQVsAGNgZ15OUBMzEU6ekohy4b1hC9fTtq1rYnv3gfy+LNld3Y5ki5xqlvEnE43ZnMGATbj\nMmTnsNst5b4tsjNERBr3e1GS2iHAloI9JDBPbeGo0ZC9PQbQcMX7VH9nIf8MGFqgDC5nKqdOOXC7\nS84gS18+PPuafhZwtVLqB2A98AwwHHheKbWFzE7i5VrrE2ReKtqctd54rbXbx2OWTqmphIwfi9ds\nxvHKTLknQAgD7bpvJCnhlWjywXxsJ48aHadE8OlMQGvtBPpd4KUbL7DuEmCJL8cpC2wL5uK/LxrX\nQw+T0ay50XGEKNfSQkLZ8dAYWv3fOK5c8Apbn5ttdCTDycfSImQ+fIjgmdPxRETievpZo+MIIYAD\nnW/nVNOrqLH5Gyr/tsXoOIaTIlCEQp4bjyk5maTnXsAbFm50HCEEgNnM9hET8ZrNXPXGVExp5esK\n9fl87Rgu83LekegL2+YfCFi1guSW1xDT9VZITMh9o/M4HIl48fqcQQhxYQn1mxDdvT8NvvyARp+9\ng+43xOhIhpEicBEORyLrft6LLSg439ua0tO5fcIEPGYz6wc+TuyuEz5liD11gsjKlQmwFWwomxDi\nv3beN5KaG9fQ9P0FHOjck5RKlY2OZAgpApdgCwomKDj/Q64afP4u4Yf3E31bP1KaX0OQj8d3OZN8\n3FIIkZs0exh/DR7FNbMn0XzJTLY99bLRkQwhfQKFzBofy+XvzsUdEsrOwY8bHUcIcQn7u/YmrkET\n6qxfQcXdfxodxxBSBApZs3fmYHU62HXPY7jDKhgdRwhxKX5+/DFsPABXznupXE43LUWgEIVF76be\nV/8joXZ9onv0NzqOECIPTjW/hoM3dKOS3kHt9V8aHafYSREoLF4vV817CZPXyx+PjMPrbzE6kRAi\nj3YMGUt6QCDN33oVf5fT6DjFSopAIblsyzoi//qVI206cfJqmSZaiNIkuXI19F0PYYs9hfp4sdFx\nipUUgUJgdru5YvH/4fHzZ8eQMUbHEUL4QPe5H1dEFdSnS8vVvEJSBApBgxXLCDl2iL23DyCpRt6m\niRZClCwZtiB23j8KP3cqVyyZaXScYiNFoICs8bE0fX8+bnsYfw981Og4QogCOHDT7cQ2akat71ZT\ncfcfRscpFlIECujy9+ZicSWxa9CjpIXK/EBClGpmM38+/DQALRa+At6yP22LFIECsB+Mpt7q/+Go\nUYfoHgOMjiOEKASnml/DoetvIeLvP6ixcY3RcYqcFIECuOLNVzF7Mtjx0BgZEipEGfLXQ6Px+Fto\n/tZrmN1le5ZRKQI+itjxC9V/+o6YZldztE0no+MIIQqRs1pN9vYYQMjxw9Rf+YHRcYqUFAFfeDy0\nWDQDgD+HPgUmk8GBhBCFbffAR3AH22n6wQIsjvxPBV9aSBHwQc2Na6i4ZycHb+hGXOMrjI4jhCgC\n7tAK/DNgKFZHAo0/WmR0nCIjRSCfzG43zd96DY+/hZ0PPGF0HCFEEYrqdQ/OytVo+MUygk4cMTpO\nkfD5eQJKqWeA2wELMA/4AVgKeICdWuvhWesNAYYCacBUrfXqAmY2VP2VHxB84gh77rgXZ7WaRscR\nQhQhjzWAnYNH0Xr60zR7eza/PDPd6EiFzqczAaXUDUAbrXVb4EagFjATGK+1vgEwK6V6KqWqACOA\nNkBXYJpSqtQOo/F3OmjywQLcwXZ2D3zE6DhCiGJwsFN34uo3odZ3qwiL/sfoOIXO18tBtwA7lVJf\nAF8Cq4CWWutNWa+vAboArYDNWut0rXUiEAWU2ovo6n9LCHAkoO96CHeoPCtAiHLBbOavB5/E5PXS\n/K2yN52Er0UgArga6AMMA94/b18OIBSwAzm71ZOAMB+PaajA0ydp9Nm7JFeMJOqOe4yOI4QoRieu\nbsfJFq2ptm0TETt+MTpOofK1T+A0sFtrnQ7sUUqlADVyvG4H4oFEMovB+ctzFRmZ/2f7Fiar1UNI\ncCzBIYEANJ2/GP/UZPTjEwiKKJ6zgGSnFQB7VgYjJDutmM0WQzNk5wBpi+wMYFxblKR2yM5QHFmi\nRzxD5Yfu4Kq3X+OnRZ9hxk1EhJ2wMGPfqwrK1yKwGRgJvKaUqg4EA98qpW7QWm8EugEbgG3AVKWU\nFbABjYGdeTlATIzDx2iFIzHRQZIzFQ8pBB85QI0vP8JxWW1239gDb1JKsWRwOt3Y7RYcxXS8i2Uw\nmzMIsBmXITuHtMXZDBGRGNYWJakdAmwp2EMCi6UtHLUac7j9zdTY/A32tas4cdV1nDrlwO0uOYMs\nffnw7FP6rBE+vyulfgFWkHlJaDTwvFJqC5kjhpZrrU8Ac8gsGuvJ7DgudfdgN1s6G3NGOjvvHyXT\nQwhRjv11/yg8Zj+avz0LU0a60XEKhc9DRLXWz1xg8Y0XWG8JsMTX4xgtfO/f1Nq4hthGzTh8/S1G\nxxFCGCipZl3+veVO6q35hPrfr4ErhxodqcBKznlMCXX5O3OAzE8AMj2EEOLvQY+SYbHS4uPFkJpq\ndJwCkyJwCZH/7KD6zxs5ecW1nGzZ1ug4QogSIDmyKtE9BhASc5zQjz80Ok6BSRG4GK+XlsvmAbBT\nzgKEEDns7j+UtMAgKsybA06n0XEKRIrARdi2bKLqru0ca9WB05e3NDqOEKIEcYdX5O8e/fE/FYNt\nSemeXE6KwIV4vVScmTlHyM7BjxscRghREu3qOZCM0DCC5r6GKSFPtz+VSFIELsD69VcE7viTf9ve\nRHyDpkbHEUKUQGnBduKHDsMcH49t/lyj4/hMisD5PB6Cp7+E12zmj/6lf/iXEKLoJNx7P56ISGyL\n5mOKPW10HJ9IETiPdfVK/Hf9RVKPXiTUrGt0HCFECeYNCsI18gnMSQ6CSunZgBSBnDwegme8hNfP\nj7gRo4xOI4QoBZLve5CMylWwLV6A6dQpo+PkmxSBHAJWfIb/P7tJ7duftDpyFiCEyAObDdeo0Zhc\nToLemG10mnyTIpAtI4OgGdPw+vnhfPIpo9MIIUqRlEGDyahWHdtbizCdPGl0nHyRIpAl4LNP8N8b\nRcqAQXjkLEAIkR+BgbhGjcGUnEzQ668ZnSZfpAgApKcT9OoreC0WXE+MNTqNEKIUShl4Lxk1amJ7\nZwnmE8eNjpNnUgSAgM+X478vmpT+g/DUrGV0HCFEaWS14np8NKaUFGxzS0/fgBSBjAyCXpuB198f\n1+NPGp1GCFGKpQwYRMZlNbC9+1ap6Rso90Ug4ItPM/sC+g/EU6u20XGEEKWZ1Ypr5JOZfQOlZKRQ\n+S4CGRkEzZyedRYw2ug0QogyIOXue8iofhm2pW9iiokxOk6uynURCPjyc/yj9pDS7248tesYHUcI\nURYEBJw9G5g3x+g0uSq/RcDjyTwL8POTswAhRKFKGXhv5n0Dby8u8XcRl9siYF39Jf76H1L79pf7\nAoQQhSsgANfIJzC5XAQtfMPoNJfk84PmAZRSlYFfgc5ABrAU8AA7tdbDs9YZAgwF0oCpWuvVBTlm\nofB6CZ45A6/ZjGuUnAUIIQpfyt33EvTa/xG4ZBGu4SPxhlcwOtIF+XwmoJTyBxYArqxFM4HxWusb\nALNSqqdSqgowAmgDdAWmKaUsBcxcYNZ1X+O/6y9Se/Umo14Do+MIIcoim43kR0diTnJge3Oh0Wku\nqiCXg/4PmA8cBUxAS631pqzX1gBdgFbAZq11utY6EYgCrijAMQvO6yXotRkAuEaNMTSKEKJsS773\nfjwVK2JbNA9TksPoOBfkUxFQSg0GTmqt15FZAM7flwMIBexAQo7lSUCYL8csLJYfvsfy26+k3nY7\nGY2bGBlFCFHWhYSQ/PBwzPHxBL69xOg0F+Rrn8D9gEcp1QVoAbwLROZ43Q7EA4lkFoPzl+cqMtLu\nY7RczJ0JQMALky55DKvVQ0hwLMEhgUWTIw+SnVYA7AZnMJsthmbIzgHSFtkZwLi2KEntkJ3BiCxm\n3ERE2AkLy+W96unRMG8OIQvnEvLMaAgKKp6AeeRTEci67g+AUmoD8AgwQynVQWv9A9AN2ABsA6Yq\npayADWgM7MzLMWJiCv/Uyf+nrVTYuJHUzjeTWLMhXOIYiYkOkpypeEgp9Bx55XS6sdstOJKMzWA2\nZxBgMy5Ddg5pi7MZIiIxrC1KUjsE2FKwhwQa0hYuZyqnTjlwu3O7oGIm6MGhBM+cQdJrr5M89NEi\ny+TLh+fCHCI6BnhBKbUFsADLtdYngDnAZmA9mR3H7kI8Zr4Ez8ruC5CZQoUQxSd5yKN4g4KwzXsd\n3Ia9BV5QgYaIAmitO+X49sYLvL4EMPximP9ff2LdsB53m3akt2ptdBwhRDnirVSJ5HvuJ2jhGwQu\n/5iUu+8xOtIZ5eZmMduczAc9yEyhQggjJA97DK/Fgu311yAjw+g4Z5SLIuC3by8BK78grXkL0jp2\nNjqOEKIc8lS/jJS7BuAfvRfrVyuNjnNGuSgCtrmzMXk8JI98Akym3DcQQogikPzY43hNJoJmzwSv\n1+g4QDkoAuZjRwn8+APS69UntXtPo+MIIcqxjPoNSe3RC8uOP7B8v8HoOEA5KAK2+XMxpaWR/Ngo\n8PMzOo4QopxLzuqXDJoz0+Akmcp0ETDFxxH43lIyqlYjpW9/o+MIIQTpzVvg7ngT1i2b8P9tm9Fx\nynYRsL39JmZnEskPD4eAAKPjCCEEAK4RTwAQVAIeSF92i0ByMrbFC/CEhpFy72Cj0wghxBlp7a4n\n7aqWWL9aiV90lKFZymwRCPzfh5hPxZAy+EG89tDcNxBCiOJiMuF6bBQmrxfbvLmGRinwHcNFZe++\nf0lP9/GGiowMrpn9Kh6LhR2dupK2Jzrfu0hKcuByeQkKLqKJ7IQQ5Zr71h6k161H4P8+wPnUeLxV\nqhiSo8QWgajD8dhCI3Nf8QIu+2EttsOH2NetL4eCap597E0+uJLTiE+IISKysk8ZhBDikvz8SH50\nJPaxowh6cwHOZycZEqPsXQ7yemn8yRK8JhO67/1GpxFCiItKuWsAnohIAt9+07CHzpTYMwFfRe7Y\nRkX9F4fbdyGphjxAXghRNLxeLw5HYoH3Y773firNnA6L5pP40MP53t5uD8VUgJkQylwRaLT8LQB0\nHzkLEEIUnWSXk43bYwmvWKlA+7FedTN9Al8n8M3FfHX1rXj98/62nOxy0qV1A0JDfX9gY5kqAvaD\n0VT/eSOnml5FbNOrjI4jhCjjAm1BBR88Emzn31t603DFMhr9toVDnboXTrg8KlN9Ao0+XQrIWYAQ\nonSJuvNevGYzavlbxT6xXJkpAgGxMdRevwJH9VocbdMp9w2EEKKEcFaryeH2N1Nh724i//i5WI9d\nZopAgy8/wC8tjT29B8tEcUKIUmdP78EAqE/fLtbjloki4JeSTP2VH5IaGs6BLr2MjiOEEPkW26QF\nMc2uptovPxD6b/FNJVEmikCdbz4nwJFAdI8BZATajI4jhBA+2ZPVn5ndv1kcfBodpJTyB94C6gBW\nYCrwN7AU8AA7tdbDs9YdAgwF0oCpWuvVBU6dk8dDw8/eIcNiZe/tdxfqroUQojgdva4jjstqU2vD\nSv564AlSK0QU+TF9PRMYBJzSWncAugJzgZnAeK31DYBZKdVTKVUFGAG0yVpvmlLKUgi5z6j28/fY\njx7kYKfuxdJgQghRZMxmou64F7+0NOqv/Kh4Dunjdv8DJmZ97QekAy211puylq0BugCtgM1a63St\ndSIQBVxRgLz/0eizdwDYc+d9hblbIYQwxL9deuG2h1F/5YeY3alFfjyfioDW2qW1diql7MAnwLNA\nzvuWHUAoYAcScixPAny/te084Xv/pvKfv3DiqjYk1m1UWLsVQgjDZNiC2HdrXwITYqn17coiP57P\ndwwrpWoCnwFztdYfKaWm53jZDsQDiWQWg/OX5yokJJCgkMBLrtP0y2UAHBo4BHsu6+aXGTe2IGuh\n7zc/kp1WAMMzmM0WQzNk5wBpi+wMYFxblKR2yM5gRJaibIdjAx6g0fKlNP7iPU73GQgXmRvIjJuI\nCDthYb7ftexrx3AVYC0wXGv9Xdbi35VSHbTWPwDdgA3ANmCqUsoK2IDGwM68HCMpKYUMc8pFXw88\nfZJq61eRWLMe+5u1hqSLr+sLlzOVZJcbRyHvNz+cTjd2u8XwDGZzBgE24zJk55C2OJshIhLD2qIk\ntUOALQV7SKAhbVGU7eAIqsChDl2p/d0qbD9s4OTV7S64nsuZyqlTDtzuzIs6kZH5Lwa+9gmMA8KB\niUqp75RSG4AJwAtKqS2ABViutT4BzAE2A+vJ7Dh2+3jMc9T/8gPM6WlE3XEvmMvESFchhDgjqndm\nP2d2v2dR8elMQGs9Chh1gZduvMC6S4AlvhznYsypKdRf/TGp9jAOdL69MHcthBAlQlyjZpk3j23b\nhP1gNI5a9YvkOKXyI3TtDSsJSIxn32395OYwIUSZFXXnvQA0+GJZkR2j9BUBr5eGn7+Hx+xHdI8B\nRqcRQogic7RNJ5yVq1Fn3QosjoTcN/BBqSsCkX/+Qti/URy+/maSI6saHUcIIYqM18+fvbcPxD81\nmbprPyuSY5S6ItDw8/cAiLrjHoOTCCFE0dvfrQ/pATYarHgfMjIKff+lqggEHztE9Z82EKuaE9vk\nSqPjCCFEkUuzh3Ggy+0EnzhC9a0bCn3/paoINFjxPiavl6iegy5684QQQpQ1UT0HAdDwi/cKfd+l\npgj4JTup+/WnJFeM4NANXY2OI4QQxcZRuwHHW7al8o5thEX/U6j7LjVFoM66FVhcSUR374/XYjU6\njhBCFKvsftDCPhsoHUXA66XBlx/g8bew77Z+RqcRQohid/zaDiRVr0Wt71ZjTYwrtP2WiiJQ+fet\nhB6M5lCHrvLMACFE+WQ2s7fHAPzcqdT9+tPC222h7akINVjxAQB7ew40OIkQQhjn31vuJD3ARv2V\nHxbacNESXwSCjh+h+s/fEduoGbGNC/V5NEIIUaqkhYRyoPPtBJ84SvWfvy+UfZb4IlB/1YeYPJ7M\nswAZFiqEKOeyn6XeYMX7hbK/El0EzKkp1F2znNSwChy6oZvRcYQQwnCJdRtxskUrqvy+lbBD+wu8\nvxJdBGptWEWAI4F93frisQYYHUcIIUqEvbdn9o82/uqTAu+r5BaBrGGhXrOZ6O79jU4jhBAlxtG2\nnXBFVqP+96sxORwF2leJLQIRe3ZRIXo3R9p0IrlyNaPjCCFEieH18yf6truwpCRj/6Jgw0VLbBFQ\na78AIDqrE0QIIcRZ+7v1IcPfn7D33wOv1+f9lMwiEBNDnR+/I7FGXU5eeZ3RaYQQosRJrRDBgTad\nsO7dg2XrFp/349MzhvNDKWUC5gEtgBTgIa31vktu9NZb+KWnZT45TIaFCiHEBemufai36RsC336T\ntLbtfdpHcZwJ9AICtNZtgXHAzFy3WLCAtIBADnTpWdTZhBCi1DrZpAWpqjEBq7/EfOK4T/sojiLQ\nHvgaQGv9M3BNrlv8+y/7r+9CWkhoEUcTQohSzGQicdB9mNLTCXxvqU+7KI4iEArkfEJyulIq1+Pq\nW3oVXSIhhCgjHLffgSfETuC7b/u0fZH3CQCJgD3H92atteeSW7RtS0y1avglnCzSYJficibhTnXi\nchZsDG5BpCQ78feHDI9x/SIpyU7MZn9D2yE7h7TF2QzOpERczlRDM5SEdnA5HZhxG9IWJaEdkl1O\nvMHVSO03ANuSRT7tw+QtwNCivFBK3Ql011o/oJS6Dpiotb6tSA8qhBAiT4rjTOBzoItSKnsM0/3F\ncEwhhBB5UORnAkIIIUquknmzmBBCiGIhRUAIIcoxKQJCCFGOSREQQohyrDhGB11QbnMKKaV6ABOB\nNOBtrfWbhgQtBnloiwHA42S2xV9a60cNCVoM8jrXlFJqIXBaaz2+mCMWmzz8XlwLvJr17XFgkNba\nXexBi0Ee2mIg8CSQTub7xQJDghYjpVRr4GWtdcfzlufrvdPIM4GLzimklPLP+r4zcCMwVCkVaUTI\nYnKptggEXgBu0FpfD4QrpbobE7NY5DrXlFLqYaBZcQczQG5tsQgYrLXuQObULLWLOV9xyq0tZgCd\nyJymZrRSKqyY8xUrpdRYYDEQcN7yfL93GlkELjWnUBMgSmudqLVOAzYDHYo/YrG5VFukAm211tm3\nRPqT+UmorLrkXFNKqTbAtcDC4o9W7C7aFkqpRsBp4Eml1PdARa11lBEhi0luc5D9CVQAbFnfl/Wx\n73uBOy6wPN/vnUYWgUvNKXT+aw6gLFf2i7aF1tqrtY4BUEqNAIK11usNyFhcLtoWSqmqwCTgMaA8\nzDF+qb+RCKANMIfMT32dlVI3Fm+8YpXbHGS7gN+Av4BVWuvE4gxX3LTWn5N56et8+X7vNLIIXGpO\noUQyf5hsdiC+uIIZ4JLzKymlTEqpGcBNwJ3FHa6YXaot+gKVgK+AZ/6/nTvWhSiIwjj+19AodEoi\nklPRiEYiEaVtJDoFIRG1t6AULY2C2jMICqVCPsEDSBQaiUSyirnFjWTvWsUsO9+v250tTk7uzrlz\n5t4B1iNiI3N8OTXl4hV4lPQg6ZN0l9z9hN7/q2MuImIGaJHaYZPAeESsZY/wb+h57uxnEbgCVgCq\nM4XuamP3wHREjEXEMGk5c5M/xGyacgGp9zsiabXWFhpUHXMh6UjSvKRlYB84k3TanzCzaLounoHR\niGilS9AAAAC3SURBVJiqPi+S7oYHVVMu3oB34ENSG3ghtYZK8H1F3PPc2bdjI2q7/bPVV1vAHKnd\ncRwRLdLSfwg4GeTd/qZckJa4t8BlNdYGDiVd5I4zh27XRe13m0AU8nRQp//IEnBQjV1L2ssfZR4/\nyMUusE3aQ3sCdqoV0sCKiAngXNJC9QThr+ZOnx1kZlYwvyxmZlYwFwEzs4K5CJiZFcxFwMysYC4C\nZmYFcxEwMyuYi4CZWcFcBMzMCvYF/cT8AtnFNQAAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline \n", "import numpy as np\n", "import matplotlib.pylab as plt\n", "import pandas as pd\n", "\n", "pd.set_option('display.width', 500)\n", "pd.set_option('display.max_columns', 100)\n", "import seaborn as sns\n", "\n", "\n", "\n", "\n", "## FUNCTIONS \n", "# target distribution p(x) \n", "p = lambda x: 6*x*(1-x)\n", "\n", "# number of samples\n", "n = 10000\n", "\n", "sig =100\n", "\n", "\n", "#intitialize the sampling. Start somewhere from 0..1\n", "x0 = np.random.uniform()\n", "\n", "\n", "x_prev = x0\n", "\n", "x=[]\n", "k=1\n", "i=0\n", "while i 1): # MAKE SURE WE STAY WITHIN BOUNDS\n", " x_star = np.random.normal(x_prev, sig)\n", " \n", " \n", " \n", " P_star = 6*x_star*(1-x_star) #p(x_star);\n", " P_prev = 6*x_prev*(1-x_prev) #p(x_prev);\n", " U = np.random.uniform()\n", " \n", " A = P_star/P_prev\n", " if U < A:\n", " x.append(x_star)\n", " i = i + 1\n", " x_prev = x_star\n", " else :\n", " x.append(x_prev)\n", " x_prev = x[i] \n", " i = i + 1\n", " \n", " k=k+1\n", "\n", "\n", "\n", "e,q,h=plt.hist(x,10, alpha=0.4, label=u'MCMC distribution') \n", "\n", "\n", "\n", "xx= np.linspace(0,1,100)\n", "plt.plot(xx, 0.67*np.max(e)*p(xx), 'r', label=u'True dsitribution') \n", "plt.legend()" ] } ], "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.12" } }, "nbformat": 4, "nbformat_minor": 0 }