09_02_random_distributions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Non-uniform random distributions\n",
    "--------------------------------\n",
    "\n",
    "In the previous section we learned how to generate random deviates with\n",
    "a uniform probability distribution in an interval $[a,b]$. This\n",
    "distributioon is normalized, so that $$\\int _a^b {P(x)dx}=1.$$ Hence,\n",
    "$P(x)=1/(b-a)$.\n",
    "\n",
    "Now, suppose that we generate a sequence $\\{x_i\\}$ and we take some\n",
    "function of it to generate $\\{y(x_i)\\}=\\{y_i\\}$. This new sequence is\n",
    "going to be distributed according to some probability density $P(y)$,\n",
    "such that $$P(y)dy=P(x)dx$$ or $$P(y)=P(x)\\frac{dx}{dy}.$$\n",
    "\n",
    "If we want to generate a desired normalized distribution $P(y)$, we need\n",
    "to solve the differential equation: $$\\frac{dx}{dy}=P(y).$$ But the\n",
    "solution of this is $$x=\\int _0^y {P(y')dy'}=F(y).$$ Therefore,\n",
    "$$y(x)=F^{-1}(x),\n",
    "$$ where $F^{-1}$ is the inverse of $F$.\n",
    "\n",
    "### Exponential distribution\n",
    "\n",
    "As an example, let us take $y(x)=-\\ln{(x)}$ with $P(x)$ representing a\n",
    "uniform distribution in the interval $[0,1]$. Then\n",
    "$$P(y)=\\frac{dx}{dy}=e^{-y},$$ which is distributed exponentially. This\n",
    "distribution occurs frequently in real problems such as the radioactive\n",
    "decay of nuclei. You can also see that the quantity $y/\\lambda$ has the\n",
    "distribution $\\lambda\n",
    "e^{-\\lambda y}$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHtBJREFUeJzt3Xt4FfW97/H3F4Lhfg8IBA3KTbwUNV5qrReoCmrFbrFibUs9nIN2l11t3bubarda2nrU7l1bW+upj7ZaWxUOxYpCtYJW1K1AAOUWhIAoIGq4X0KAhO/+YwbWSkzIgqxkkpnP63nWk/n95rdWvsPwfNbKb2bNmLsjIiLJ0CLqAkREpPEo9EVEEkShLyKSIAp9EZEEUeiLiCSIQl9EJEEU+iIiCaLQFxFJEIW+iEiC5ERdQHXdu3f3goKCqMsQEWlWFixYsMnd8+oa1+RCv6CggKKioqjLEBFpVszsg0zGaXpHRCRBFPoiIgmi0BcRSRCFvohIgij0RUQSRKEvIpIgCn0RkQRR6IuIJIhCX0QkQWIX+gUTZ0RdgohIkxW70BcRkdop9EVEEkShLyKSIAp9EZEEySj0zWyEmb1nZiVmNrGG9ReY2UIzqzCz0dXWjTWzVeFjbLYKFxGRI1dn6JtZS+AhYCQwBLjezIZUG/Yh8C3gqWrP7QrcBZwDnA3cZWZd6l92LTZvZlDp2gZ7eRGR5i6TT/pnAyXuvsbd9wHPAKPSB7j7WndfDByo9tzLgJfdfYu7bwVeBkZkoe6qPv4YevaE7t350zM/yvrLi4jERSah3wdYl9ZeH/ZlIqPnmtl4Mysys6LS0tIMXzpNjx6wcycAeWXbYPPmI38NEZEEaBIHct39EXcvdPfCvLw6b/H4WS1awODBqXZxcfaKExGJkUxCfwPQN62dH/Zloj7PPTInnZRaVuiLiNQok9CfDwwws35mdgwwBpie4eu/BFxqZl3CA7iXhn3ZNyTt2PLy5Q3yK0REmrs6Q9/dK4AJBGFdDExx92VmNsnMrgIws7PMbD1wLfA7M1sWPncL8BOCN475wKSwL/v0SV9EpE7m7lHXUEVhYaEXFRUd+ROLi1Of9vv2hQ8/zG5hIiJNmJktcPfCusY1iQO5WdG/P+TkBMvr1h06m0dERFLiE/qtWsGAAYeaV014NMJiRESapviEPlQ5mNt/87rDDBQRSaZ4hX7awdwBmxT6IiLVxSv09UlfROSw4hX6aZ/0+2/W2TsiItXFK/QHDeIABsBx2z6B8vKICxIRaVriFfpt2rCuc08AWvoBWLky4oJERJqWeIU+sKpb2qV+9M1cEZEqYhf6Jemhr2vwiIhUEb/Q735cqqFP+iIiVcQv9PVJX0SkVvEO/ZUroaIiumJERJqY2IX+rty2fNShe9DYvx9WrYq2IBGRJiR2oQ/wXt7xqcaSJdEVIiLSxMQy9Ffk9Us1Fi+OrhARkSYmlqFf3KMg1VDoi4gcEsvQX5FXkGoo9EVEDoll6K/pms++FuFdtD74ALZvj7YgEZEmIpahX9Eyh9Xd8lMdOpgrIgLENPQBinukHcxV6IuIADEOfc3ri4h8lkJfRCRBYhv6n5neOXAgumJERJqI2IZ+absubGnTMWjs3BmcxSMiknCxDX3MNMUjIlJNfEMfzeuLiFQX69CvcjkGnbYpIhLv0NeF10REqop16K/q3pdKCzdx1SooK4u2IBGRiGUU+mY2wszeM7MSM5tYw/pcM5scrp9rZgVhfysze8LMlphZsZn9MLvlH155q9as7dI7aBw4AEuXNuavFxFpcuoMfTNrCTwEjASGANeb2ZBqw8YBW929P/AAcF/Yfy2Q6+6nAmcCNx18Q2gsy3qekGosWNCYv1pEpMnJ5JP+2UCJu69x933AM8CoamNGAU+Ey1OB4WZmgAPtzCwHaAPsA3ZkpfIMLenZP9VYuLAxf7WISJOTSej3AdaltdeHfTWOcfcKYDvQjeANYDewEfgQ+E9331LPmo/I0mNPTDX0SV9EEq6hD+SeDVQCvYF+wG1mdkL1QWY23syKzKyotLQ0qwUs65kW+kuXwt69WX19EZHmJJPQ3wD0TWvnh301jgmncjoBm4GvAS+6+353/xR4Eyis/gvc/RF3L3T3wry8vCPfisPY0bo9H3Q+Nmjs36+DuSKSaJmE/nxggJn1M7NjgDHA9GpjpgNjw+XRwCvu7gRTOsMAzKwdcC6wIhuFH4kq8/qa4hGRBKsz9MM5+gnAS0AxMMXdl5nZJDO7Khz2GNDNzEqA7wMHT+t8CGhvZssI3jz+4O6N/i2ppcemQv/PDz/b2L9eRKTJyMlkkLvPBGZW67szbbmc4PTM6s/bVVN/Y1uaNq9/yserI6xERCRasf5G7kHpZ/AMLn0f9u2LsBoRkegkIvS3tenIuk49AcitrIDlyyOuSEQkGokIfag6xaODuSKSVIkJ/SXH6gweEZHEhH6VL2npcgwiklCJCf0qn/TffRcqKqIrRkQkIokJ/S1tO7GhQ/ht3/JyWLYs2oJERCKQmNAHWNxrQKoxb150hYiIRCRRob+o96BU4+23oytERCQiCn0RkQRJVOgvObY/FQfvmbt8OWzbFm1BIiKNLFGhX96qNcU9+qU65s+PrhgRkQgkKvQBFvUenGpoikdEEiaBoa95fRFJLoW+e3TFiIg0ssSF/touvaFr16CxZQuUlERbkIhII0pc6GPGK53T7s2uKR4RSZDkhT6a1xeR5Epo6OsMHhFJpkSG/ru9B4JZ2HgXysqiLUhEpJEkMvR35raDk04KGpWV+pKWiCRGIkMfgPPOSy2/8UZ0dYiINKLkhv4FF6SW58yJrg4RkUaU3ND/4hdTy2++qTtpiUgiJDf0jz8e+vYNlnfvhkWLoq1HRKQRJDf0zTTFIyKJk9zQB4W+iCROYkO/YOIMhs8/kOp4/XU4cKD2J4iIxEBiQx9gddd8NrXtFDS2boVly6ItSESkgSU69DFjXv7Jqfbrr0dXi4hII0h26APz+p6SamheX0RiLqPQN7MRZvaemZWY2cQa1uea2eRw/VwzK0hbd5qZvWVmy8xsiZm1zl759Te/b9on/TlzdFMVEYm1OkPfzFoCDwEjgSHA9WY2pNqwccBWd+8PPADcFz43B/gTcLO7nwxcBOzPWvVZUJxXAB07Bo2NG2H16kjrERFpSJl80j8bKHH3Ne6+D3gGGFVtzCjgiXB5KjDczAy4FFjs7u8CuPtmd6/MTunZcaBFSzj//FTHq69GV4yISAPLJPT7AOvS2uvDvhrHuHsFsB3oBgwE3MxeMrOFZvaD+pfcAIYNSy3PmhVdHSIiDSynEV7/fOAsoAyYbWYL3H12+iAzGw+MBzjuuOMauKQaXHJJann27OB8/RaJP8YtIjGUSbJtAPqmtfPDvhrHhPP4nYDNBH8VzHH3Te5eBswEzqj+C9z9EXcvdPfCvLy8I9+K+jrlFOjRI1jevBneeafxaxARaQSZhP58YICZ9TOzY4AxwPRqY6YDY8Pl0cAr7u7AS8CpZtY2fDO4EFiendKzqEULGD481dYUj4jEVJ2hH87RTyAI8GJgirsvM7NJZnZVOOwxoJuZlQDfByaGz90K/ILgjeMdYKG7z8j+ZmRB+hSPQl9EYiqjOX13n0kwNZPed2facjlwbS3P/RPBaZtN25e+dGix/NXXaF1eDq2b1FcKRETqTUcrD+rbFwYNAqB1xb7gxioiIjGj0Ce44mbBxBlVPu3/9kf/L8KKREQahkI/zf/Z2OXQ8vlrdQaPiMSPQj/N28edSoUF/ySnfLw6OH1TRCRGFPppdua2491eAwFogQdf1BIRiRGFfjVvFJyeasycWftAEZFmSKFfzez+Z6UaM2fqFooiEisK/WqWHNuf0nadg0ZpKcyfH21BIiJZpNCvxq0Fr55QmOp44YXoihERyTKFfg1eOTFtimdG07xqhIjI0VDo1+CNgtPZ1yK8QsWiRbCh+kVFRUSaJ4V+DXbltmVu+g3TdRaPiMSEQr8Wr6ZP8WheX0RiQqFfiyqnbs6aBeXl0RUjIpIlCv1afNClNwwMvp1LWZlumC4isaDQP5wrr0wt//Wv0dUhIpIlCv3D+ad/Si0/+yxUVkZXi4hIFij0D6Pf9C180r5r0Cgthddfj7YgEZF6UugfhlsLXhz4+VTH1KnRFSMikgUK/Tq8OPALqca0aboAm4g0awr9OszrezLk5QWNjRvhrbeiLUhEpB4U+nWobNESrr461fGXv0RXjIhIPSn0MzF6dGp56lRwj64WEZF6UOhn4uKLoUt40/R163SNfRFpthT6mWjVCkaNSrUnT46uFhGRelDoZ2js7n6pxtNP64taItIsKfQz9EbBUErbhrdR3LgRXnkl2oJERI6CQj9DlS1a8vxJF6Q6/vSn6IoRETlKCv0jMO2UYWmNabB7d3TFiIgcBYX+EVja80QYPDho7NoFzz0XbUEiIkdIoZ+BgonhzdHN4OtfT63QFI+INDMZhb6ZjTCz98ysxMwm1rA+18wmh+vnmllBtfXHmdkuM/vX7JQdnfPX90o1/v53+OST6IoRETlCdYa+mbUEHgJGAkOA681sSLVh44Ct7t4feAC4r9r6XwB/q3+50VvfqSdz808OGpWVTPraf0RbkIjIEcjkk/7ZQIm7r3H3fcAzwKhqY0YBT4TLU4HhZmYAZnY18D6wLDslR++vJ198aPmri/8O7hRMnJGaBhIRaaIyCf0+wLq09vqwr8Yx7l4BbAe6mVl74N+BHx/uF5jZeDMrMrOi0tLSTGuPzPMnXUBZq1wABm/6AN5+O+KKREQy09AHcu8GHnD3XYcb5O6PuHuhuxfmHbyMcRO2K7ct00+6MNXxyCPRFSMicgQyCf0NQN+0dn7YV+MYM8sBOgGbgXOA+81sLXArcLuZTahnzU3C05+7LNWYPJmO5Yd9XxMRaRIyCf35wAAz62dmxwBjgOnVxkwHxobLo4FXPPBFdy9w9wLgl8A97v6bLNUeqXd7DWR5j/B6PHv2MGr5PyKtR0QkE3WGfjhHPwF4CSgGprj7MjObZGZXhcMeI5jDLwG+D3zmtM7YMeOpoSMPNb/2zou6zr6INHk5mQxy95nAzGp9d6YtlwPX1vEadx9FfU3ac0Mu5PZXH6Pt/r2cVLqWoRtXRl2SiMhh6Ru59bAztx3PD05dhO0bi3TKpog0bQr9enpq6IhDy1cWz4GPP46wGhGRw1Po19O7vQfB5z8PQG5lBfz2txFXJCJSO4V+Nnzve6nlhx+GPXuiq0VE5DAU+tnwla+wvmOPYHnTJnjqqWjrERGphUI/Cwp+9BKPn3llquOBB3T6pog0SQr9LJly2qXsbtU6aCxbBrNmRVuQiEgNFPpZsqN1e6acdkmq4+c/j64YEZFaKPSz6PeFo6i08J/05Zdh7txoCxIRqUahn0XrOh/Lc0PSrr75059GV4yISA0U+ln20Llf5QAWNF54ARYtirYgEZE0Cv0sW929LzMHfSHVoU/7ItKEKPQbwG/Ouy7VmDYtOJsHdEtFEYmcQr8BrOjRj5f7n5PqmDQpumJERNIo9BvIg+eNSTWmTIEFC6IrRkQkpNBvIEt6DeDFgZ9PdUyM/31lRKTpU+g3oJ9f8E1oEf4Tz5rF+e/rTB4RiZZCvwGt7tYXxo071P731x7H/ECEFYlI0in0G9pdd0GbNgCc+slqrljxRsQFiUiSKfQbWp8+cMsth5oT//G4rrcvIpFR6DewgokzOG3PUDa36QhA/o5P4d57I65KRJJKod8IdrRuz/0Xjk113HcfrFkTXUEiklgK/UYy5bRLeKfXgKCxd2/VWyyKiDQShX4jcWvBXV+6OXUxtunTYebMaIsSkcRR6Deid3sPYnL6jVYmTIBdu6IrSEQSR6HfyH5+4Vjo0iVovP8+3HFHtAWJSKIo9BvZlrad4Je/THX8+tfwhs7dF5HGodCPQMGyrrx6wplBwz341q7O3ReRRqDQj4IZt182gZ3HBN/UZeVKuPPOaGsSkURQ6EdkY8c87rk4dV0e/uu/YPbs6AoSkUTIKPTNbISZvWdmJWb2mWsEm1mumU0O1881s4Kw/xIzW2BmS8Kfw7JbfvP29Ocug0vCs3nc4RvfgNLSaIsSkVirM/TNrCXwEDASGAJcb2ZDqg0bB2x19/7AA8B9Yf8m4MvufiowFngyW4XHghlnDfwG5OUF7Y0b4cYbgzcAEZEGkMkn/bOBEndf4+77gGeAUdXGjAKeCJenAsPNzNx9kbt/FPYvA9qYWW42Co+L0vZd4fHHUx0zZsCDD0ZWj4jEWyah3wdYl9ZeH/bVOMbdK4DtQLdqY64BFrr73qMrNb4K5jiPFqa9j/7bv+k0ThFpEI1yINfMTiaY8rmplvXjzazIzIpKEzqnff+F32JJzxODxv79MHo0rF8fbVEiEjuZhP4GoG9aOz/sq3GMmeUAnYDNYTsfeBb4pruvrukXuPsj7l7o7oV5B+e3E2ZfTitu/sod0L170PHJJ3DNNVBeHm1hIhIrmYT+fGCAmfUzs2OAMcD0amOmExyoBRgNvOLubmadgRnARHd/M1tFx9WGTj1gyhRo2TLomDcPbr5ZB3ZFJGvqDP1wjn4C8BJQDExx92VmNsnMrgqHPQZ0M7MS4PvAwdM6JwD9gTvN7J3w0SPrWxEnF18cnLN/0BNPwKRJ0dUjIrFi3sQ+RRYWFnpRUdFRP79g4owsVhMRd34+81dcu3RWqu/RR6vcZF1EJJ2ZLXD3wrrG6Ru5TZEZt4/4DnMKTk/13XSTrr8vIvWm0G+i9rdsxbev/iHLepwQdFRWBmf0vPbaoTEFE2fE4y8bEWk0Cv0mbHduW7517d1w/PFBx549cMUV8MYbCnsROSoK/SautH1XLrr0dj5p3zXo2L0bRo7k9A0roi1MRJolhX4zsLZrH64fcw/07Bl07NrFH6f8B+d+uDjawkSk2VHoNxNruuXzpSvvYlPbTgB02LeHJ6bcxSWr3gY0vy8imVHoNyMl3Y/j+jH3HJrqya3cz8PP3lP1gm0iIoeh0G9mVuUdzzU33M/azr0AyPEDcOON3DbnScwPRFydiDR1Cv1maH3nYxn99ftZ3qPfob5/eWsyv/3rvZz0/b9EWJmINHUK/WZqU7suXPe1e/lHvzMP9Y1c+d9M/fMPYHWN17UTEVHoN2c7c9sxbvSdcOuth/pO/nQNnHEGTJumA7si8hkK/WauskVLeOABJl42gX0tcoLOHTvgmmu4a9bvGPSvzyr8ReQQhX4MFEycwTNDRzD66/ezrlPPQ/03Lnie5x+/lVM+LomwOhFpShT6MbK410Cu+Nav+PuAcw/1Ddz8Ic8+eRvcfTfs1Z0qRZJOoR8zO1q3Z/xX7uCOS/+ZslbBPehbHaiEH/8Yhg6tcsE20Je6RJJGoR9HZvz59MsZceNvmN9nSKp/xQq46CIYOxY2VL/jpYgkgUI/xj7s0ovrvvZ/uXv4eOjQIbXij3+EAQPgRz+i/d6y6AoUkUan0I+5Ay1a8njhVVBczIxBX0it2LMHfvYzXvvd/2bsgudh377oihSRRqPQT4iCX7/Dd67+ITdc91M4PXVHrm57dvDjWb+D/v3hwQehTJ/8ReJMoZ8wbxYMhaIibrnyNtZ3TLtH/bp1cMstwQ1bfvITTrt1sg7yisRQTtQFSOMruP1vcPLFvDjoC3xj4QvcPPcvdC/bHqzctAnuvJO3WrXmuSEX8eehI6ItVkSyytw96hqqKCws9KKioqN+vj6ZHrnc/Xv56pKXuWnuNPJ3fPrZAWedBTffDNddB+3aNX6BIlInM1vg7oV1jdP0jrC3VS5PnnElF41/hFuuvA1OPbXqgPnzYdy44M5dN9wAL7ygA78izZQ+6ctnuVO4YTk3LPobl7/3JrmV+z87pksXGD0arr4ahg2D1q0bv04ROSTTT/oKfTmsLmXbuWbpbK5b/DIDNq+rcUxZq1zaXj4Cvvxlzl7YinkPfbORqxQRhb5klzsnlb7PqOWvcWXxHPJ3lNY6dFW3vvz38acx9o7/BRdeCN26NWKhIsmk0JcGY36AMzas4LKVbzF89TxO3HKYSzqYweDBcO65weOcc+DkkyEnOHHs4P5ae+8VjVG6SGxlGvo6ZVOOmFsLFuQPYUH+EO4ZNo5+WzYwrGQew1fP58wNy8mtrEgb7FBcHDz+8Iegr107KCyEM87g2hUHWJFXEHwprG3bSLZHJEn0SV+yKnf/Xs7cUMx5Hy7mvA/e5YxPSqCysu4ntmgRfCv41FN58KMcvnvTyOD6QP37Q15e8BeDiNRK0zvSJKz9j4tg4ULuufMPDP3oPYZuXEnvnZuO7EU6dgzCv39/OPFEOO44yM+Hvn2Dn1276k1BEk/TO9IkFPzkH8HCOdcc6lv7L0ODc/8XL2bGU39ncOlaCrZupKUfqPlFduyAhQuDR01at069CfTuHfxl0KNHzT87dNAbhCRaRqFvZiOAXwEtgUfd/d5q63OBPwJnApuB69x9bbjuh8A4oBL4rru/lLXqpXnq0wf69KHg7VZw9ZkAtN5fzoBN6xi06QMKtn5EwZaPKNi2kYKtH9F+357Dv155OZSUBI+65OZC9+7QuTN06hT8DJcfWrSJnbntmHj9ual1HToExxratUs92raFY47Jwj+ESOOrc3rHzFoCK4FLgPXAfOB6d1+eNuafgdPc/WYzGwN8xd2vM7MhwNPA2UBvYBYw0N1rneTV9I5U4U63su3BG8HWjRy37WOO3bmJXjs3cUGb8uBCcbt2NX5dOTmH3gTe3+PsadWaIf17QZs2wRvLwccxx1Rt19Z3sD8n5/CPli2PbL1ZcLzk4EN/5cRWNqd3zgZK3H1N+MLPAKOA5WljRgF3h8tTgd+YmYX9z7j7XuB9MysJX++tTDdEEs6Mze06s7ldZxbkD6lxSPu9ZfTaUUrvnZvovnsb3cq20bVsO93LttO1bDvd0h5tKrJ0n+CKCti+HbZvp9/Bvk/fz85rN7TqbwQH3wzq6st0zME3lkx+HsnYo/2Zrdeq7d/yaNYdbv2ECXD55Yd/bj1kEvp9gPSvYq4HzqltjLtXmNl2oFvY/3a15/Y56mpFarArty2r8o5nVd7xdY5ts6+crnt20GHvbjrs3U3HvbvpsLeMjuW7gp97d1dZbrdvD233l9O6Yi9t9++lzb5y2u4vJ6e24w/NgXtwRlUmZ1VJ47vqqgZ9+SZxINfMxgPjw+YuM3uvHi/XHTjC00OaPW1zMmibk+Db3+7Ot799NNtc96ceMgv9DUDftHZ+2FfTmPVmlgN0Ijigm8lzcfdHgEcyKbguZlaUybxWnGibk0HbnAwNvc2ZXFp5PjDAzPqZ2THAGGB6tTHTgbHh8mjgFQ+OEE8HxphZrpn1AwYA87JTuoiIHKk6P+mHc/QTgJcITtn8vbsvM7NJQJG7TwceA54MD9RuIXhjIBw3heCgbwXwncOduSMiIg0rozl9d58JzKzWd2facjlwbS3P/Rnws3rUeKSyMk3UzGibk0HbnAwNus1N7jIMIiLScHS7RBGRBIlN6JvZCDN7z8xKzGxi1PU0BDPra2avmtlyM1tmZreE/V3N7GUzWxX+7BJ1rdlmZi3NbJGZvRC2+5nZ3HB/Tw5PMogVM+tsZlPNbIWZFZvZ5+O+r83se+H/7aVm9rSZtY7bvjaz35vZp2a2NK2vxv1qgQfDbV9sZmfU9/fHIvTDS0U8BIwEhgDXh5eAiJsK4DZ3HwKcC3wn3M6JwGx3HwDMDttxcwtQnNa+D3jA3fsDWwmu7xQ3vwJedPfBwOcItj+2+9rM+gDfBQrd/RSCE0fGEL99/Tgwolpfbft1JMFZjwMIvsv0cH1/eSxCn7RLRbj7PuDgpSJixd03uvvCcHknQQj0IdjWJ8JhTwBXR1NhwzCzfOAK4NGwbcAwgkt+QDy3uRNwAcGZcbj7PnffRsz3NcHJJW3C7/u0BTYSs33t7nMIznJMV9t+HQX80QNvA53NrFd9fn9cQr+mS0XE+nIPZlYAnA7MBXq6+8Zw1cdAz4jKaii/BH4AHLz2QTdgm7sfvEVXHPd3P6AU+EM4rfWombUjxvva3TcA/wl8SBD224EFxH9fQ+37NevZFpfQTxQzaw/8BbjV3Xekrwu/FBebU7LM7ErgU3dfEHUtjSwHOAN42N1PB3ZTbSonhvu6C8En234EV+Vtx2enQWKvofdrXEI/o8s9xIGZtSII/D+7+7Sw+5ODf/KFPz+Nqr4G8AXgKjNbSzBtN4xgrrtzOAUA8dzf64H17j43bE8leBOI877+EvC+u5e6+35gGsH+j/u+htr3a9azLS6hn8mlIpq9cC77MaDY3X+Rtir9Mhhjgecau7aG4u4/dPd8dy8g2K+vuPsNwKsEl/yAmG0zgLt/DKwzs0Fh13CCb7bHdl8TTOuca2Ztw//rB7c51vs6VNt+nQ58MzyL51xge9o00NFx91g8gMsJbvayGrgj6noaaBvPJ/izbzHwTvi4nGCOezawiuBGNV2jrrWBtv8i4IVw+QSC6ziVAP8fyI26vgbY3qFAUbi//wp0ifu+Bn4MrACWAk8CuXHb1wQ3ltoI7Cf4i25cbfsVMIIzE1cDSwjObKrX79c3ckVEEiQu0zsiIpIBhb6ISIIo9EVEEkShLyKSIAp9EZEEUeiLiCSIQl9EJEEU+iIiCfI/S58iika7ib8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot\n",
    "\n",
    "N = 10000\n",
    "r = np.random.random(N)\n",
    "\n",
    "xlambda = 0.1\n",
    "x = -np.log(r)/xlambda\n",
    "\n",
    "binwidth=xlambda*5\n",
    "pyplot.hist(x,bins=np.arange(0.,100., binwidth),density=True);\n",
    "pyplot.plot(np.arange(0.,100.,binwidth),xlambda*np.exp(-xlambda*np.arange(0.,100.,binwidth)),ls='-',c='red',lw=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "von Neumann rejection\n",
    "---------------------\n",
    "A simple and ingenious method for generating random points with a\n",
    "probability distribution $P(x)$ was deduced by von Neumann. Draw a plot\n",
    "with you probability distribution, and on the same graph, plot another\n",
    "curve $f(x)$ which has finite area and lies everywhere above your\n",
    "original distribution. We will call $f(x)$ the “comparison function”.\n",
    "Generate random pairs $(x_i,y_i)$ with uniform distribution inside\n",
    "$f(x)$. Whenever the point lies inside the area of the original\n",
    "probability, we accept it, otherwise, we reject it. All the accepted\n",
    "points will be uniformly distributed within the original area, and\n",
    "therefore will have the desired distribution. The fraction of points\n",
    "accepted/rejected will deppend on the ratio between the two areas. The\n",
    "closer the comparison function $f(x)$ resembles $P(x)$, the more points\n",
    "will be accepted. Ideally, for $P(x)=f(x)$, all the points will be\n",
    "accepted, and none rejected. However, in practice, this is not always\n",
    "possible, but we can try to pick $f(x)$ such that we minimize the\n",
    "fraction of rejected points.\n",
    "\n",
    "It only remains how to pick a number with probability $f(x)$. For this\n",
    "purpose, we utilize the method shown in the previous section, using a\n",
    "function whose indefinite intergral is know analitically, and is also\n",
    "analitically invertible. We then pick a random number $x$ and retrieve\n",
    "the corresponding $y(x)$ according to (\\[random\\_invert\\]). Then, we\n",
    "generate a second random number and we use the rejection criterion.\n",
    "\n",
    "An equivalent procedure consists of picking the second number between 0\n",
    "and 1 and accept or reject according to wether is it respectively less\n",
    "than or greater than the ratio $P(x)/f(x)$. Clearly, if $f(x)=P(x)$ all the points will be accepted.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptance Ratio:  0.171\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXJwlJQDaBiKwmCqIooBhQq7X9SRfcwI1bKFbahy21lfZW21q8trZqa7X2urTaxV/Vq9a6AK2lQt31uqESUEFWI2sAlU1WIYR87h/nkBmGQCZkkpOceT8fj3nw/X7nzMxnHkffc/I9m7k7IiKSHXKiLkBERJqOQl9EJIso9EVEsohCX0Qkiyj0RUSyiEJfRCSLKPRFRLKIQl9EJIso9EVEskhe1AWk6tKlixcXF0ddhohIizJr1qx17l5U13LNLvSLi4spKyuLugwRkRbFzJans5ymd0REsohCX0Qkiyj0RUSyiEJfRCSLKPRFRLJIWqFvZsPNbJGZlZvZxFqeP8PMZptZlZldnPLcODN7P3yMy1ThIiJSf3WGvpnlAncDZwH9gTFm1j9lsRXA14G/pby2E/Bz4GRgKPBzMzu04WWLiMjBSGdLfyhQ7u5L3L0SeBQYmbyAuy9z9zlAdcprvww86+4b3H0j8CwwPAN17239ehg7Fk46CQYMyPjbi4jERTonZ/UAVib1Kwi23NNR22t7pPna9LVtC489Brt3gxls3w5t2mT8Y0REWrpmsSPXzMabWZmZla1du7b+b1BQACUlQdsd3n8/swWKiMREOqG/CuiV1O8ZjqUjrde6+z3uXurupUVFdV46onbHHJNoL1p0cO8hIhJz6YT+TKCvmZWYWT4wGpia5vs/DXzJzA4Nd+B+KRzLvOTQX7iwUT5CRKSlqzP03b0KmEAQ1guAx919npndYGYjAMxsiJlVAKOAP5vZvPC1G4AbCX44ZgI3hGOZ169foq3QFxGpVVpX2XT36cD0lLHrktozCaZuanvtfcB9DagxPdrSFxGpU7PYkZsRqXP61alHj4qISHxCv0sX6Nw5aG/fDqvS3dcsIpI94hP6oHl9EZE6xCv0ddimiMgBxTf0taUvIrKPeIW+pndERA4oXqGvLX0RkQOKV+iXlFCZE556sGoVbNkSbT0iIs1MvEK/VSuWH9ot0V+8OLpaRESaoXiFPrCkU9KVmzXFIyKyl9iF/gedk64GodAXEdlL/EK/U9KVnHWsvojIXmIX+preERHZv9iF/l7TO4sXB7dQFBERIIahv7mwLXTtGnR27oSlS6MtSESkGYld6ANw3HGJ9rx50dUhItLMxDP0jz8+0Vboi4jUiH/ov/dedHWIiDQzsQz9C1/eVNNe8NwMiidOi7AaEZHmI5ahv7jLETXto9ZXkLe7KsJqRESaj1iG/taCNlS0LwIgv7qK4o2rI65IRKR5iGXow95b+/3WLo+wEhGR5iO+oV+UCP2j1yn0RUQgzqHfpXdNu59CX0QEiHHoL0qa3jla0zsiIkCMQ7+8cy8wA6B44xr49NOIKxIRiV5sQ39nqwLo0weAHBwWLIi4IhGR6MU29AGdmSsikkKhLyKSRbIn9HXhNRGRLAr9NLb0dY0eEYm7eId+375U5uQF7RUrYNOmAy8vIhJzaYW+mQ03s0VmVm5mE2t5vsDMHguff9PMisPxVmb2gJnNNbMFZnZNZsuvQ6tWe98+ce7cJv14EZHmps7QN7Nc4G7gLKA/MMbM+qcsdhmw0d37ALcDt4Tjo4ACdx8AnAR8e88PQlOZf1hJovPOO0350SIizU46W/pDgXJ3X+LulcCjwMiUZUYCD4TtycAwMzPAgUPMLA9oDVQCmzNSeZrmH3ZkovP220350SIizU46od8DWJnUrwjHal3G3auATUBngh+AbcAaYAXwW3ff0MCa62V+16MSHW3pi0iWa+wduUOB3UB3oAT4oZkdmbqQmY03szIzK1u7dm1GC5jXNenj3nsPdu3K6PuLiLQk6YT+KqBXUr9nOFbrMuFUTgdgPfBV4Cl33+XuHwOvAaWpH+Du97h7qbuXFhUV1f9bHMDmwrZUtD8s6FRWwvz5GX1/EZGWJJ3Qnwn0NbMSM8sHRgNTU5aZCowL2xcDL7i7E0zpnAlgZocApwALM1F4fey1ta8pHhHJYnWGfjhHPwF4GlgAPO7u88zsBjMbES52L9DZzMqBq4A9h3XeDbQ1s3kEPx73u/ucTH+JumhnrohIIC+dhdx9OjA9Zey6pPYOgsMzU1+3tbbxprLnDNt52pkrIgLE/Yzc0PyuKcfqu0dXjIhIhLIi9Fe3K2JjYbugs2kTLFsWaT0iIlHJitDHbO+duZrXF5EslR2hT8rOXM3ri0iWyprQ15a+iEgWhb629EVEsij0l3TuyY68/KBTUQHr1kVbkIhIBLIm9Hfn5LKwqDgxMGtWZLWIiEQla0IfYM7hfROdmTOjK0REJCJZFfrvdjs60XnrregKERGJSJaFftKW/ltv6cxcEck6WRX6Szr1YEt+66Dz0UfBDl0RkSySVaFfnZPLXM3ri0gWy6rQB5iTOsVDcDXOPVfkFBGJs6wL/Xe0M1dEsljWhf5eW/plZVBdHV0xIiJNLOtCf3W7ItYe0jHobNkCixZFW5CISBPKutDHbO8pHu3MFZEskn2hT8qZuZrXF5EskpWhrzNzRSRbZWXo77Uz9913ya/aFV0xIiJNKCtD/5PW7VnWsVvQqazk2I+XRFuQiEgTycrQB3ine2KKZ/DqhRFWIiLSdLI29Mt69K9pl1bMj7ASEZGmk7WhP6vnsTXt0lULdMVNEckKWRv6i7ocweb8NgB03bqBnps+irgiEZHGl7WhX52Ty9s9jqnpD9EUj4hkgawNfYCyHslTPAp9EYm/vKgLiFJZz3135iZfYnnZzec0eU0iIo0pq7f03+nWD3JzAei3bgXtd2yNuCIRkcaV1aH/aX4hnHhiTX/wqgURViMi0viyOvQBOO20mqZ25opI3Cn0k0K/VFv6IhJzaYW+mQ03s0VmVm5mE2t5vsDMHguff9PMipOeG2hmM8xsnpnNNbPCzJWfAUmhP2jNYlrt1sXXRCS+6gx9M8sF7gbOAvoDY8ysf8pilwEb3b0PcDtwS/jaPOCvwOXufhzweaB5pWr37qzo0BWAwqpKjv/wg4gLEhFpPOls6Q8Fyt19ibtXAo8CI1OWGQk8ELYnA8PMzIAvAXPc/V0Ad1/v7rszU3rmzEw6dPPkle9FWImISONKJ/R7ACuT+hXhWK3LuHsVsAnoDBwNuJk9bWazzezqhpeceW/0HljTPnXFnAgrERFpXI19clYecDowBNgOPG9ms9z9+eSFzGw8MB6gd+/ejVzSvmYckQj9IRXzaLV7F7tyWzV5HSIijS2dLf1VQK+kfs9wrNZlwnn8DsB6gr8KXnb3de6+HZgODE79AHe/x91L3b20qKio/t+igSo6dK2Z12+zaycD17zf5DWIiDSFdEJ/JtDXzErMLB8YDUxNWWYqMC5sXwy84O4OPA0MMLM24Y/B54BmeTD8jKQpns8sfzfCSkREGk+doR/O0U8gCPAFwOPuPs/MbjCzEeFi9wKdzawcuAqYGL52I3AbwQ/HO8Bsd5+W+hnNwetHJM/rz42wEhGRxpPWnL67TyeYmkkeuy6pvQMYtZ/X/pXgsM1mLXlL/6RVCyioqoywGhGRxqEzckMft+vMB516AlCwexeDV+m+uSISPwr9JDN6D6hpn6p5fRGJIYV+ktePGFTT1ry+iMRR1od+8k1T3kja0j9hzSLYti2KkkREGk3Wh36yDW06sKCoGIBW1bvhlVeiLUhEJMMU+imSj+LhmWeiK0REpBEo9FP875EnJTpPPRVdISIijUChn+KNXsezIy8/6CxYAMuXR1uQiEgGKfRT7GxVwBu9Ejt0tbUvInGi0K/F/x6ZdE04hb6IxIhCvxYvHVma6Dz3HFTqkgwiEg8K/VosPTRxC0W2boXXX4+2IBGRDFHo18Zs76N4/v3v6GoREckghf5+vKRDN0UkhhT6+zGj90DIDw/dnDMHVqXeLExEpOVR6O/H9vzW8NnPJgY0xSMiMaDQP5Czzkq0p6beIVJEpOVR6B/IyJGJ9jPPBEfyJCmeOG2vq3SKiDR3Cv0D6dMHjj8+aO/cqQuwiUiLp9Cvy/nnJ9pPPFHT1Ba+iLRECv26JIf+k0/Crl3R1SIi0kAK/boMHszqdl2C9saNurGKiLRoCv26mPFM31MS/aQpHhGRlkahn4Z9Qt89umJERBpAoZ+Gt3odDx07Bp2VK+Htt6MtSETkICn001CVmwfnnpsYmDIlumJERBpAoZ+uCy9MtB97TFM8ItIiKfTTddZZ0L590P7gAwZ++H609YiIHASFfroKC+GCC2q65y14OcJiREQOjkK/PkaPrmmeu+AVzKsBnZ0rIi2HQr8+hg2DLsGJWt22rmdIxfyICxIRqR+Ffn20agUXX1zT1RSPiLQ0aYW+mQ03s0VmVm5mE2t5vsDMHguff9PMilOe721mW83sR5kpO0JjxtQ0z174Knm7qyIsRkSkfuoMfTPLBe4GzgL6A2PMrH/KYpcBG929D3A7cEvK87cB8bj11OmnQ48eAHT+dDOfWf5uxAWJiKQvnS39oUC5uy9x90rgUWBkyjIjgQfC9mRgmJkZgJmdDywF5mWm5Ijl5MBXvlLTvWDeixEWIyJSP+mEfg9gZVK/IhyrdRl3rwI2AZ3NrC3wE+D6hpfajFxySU3zrMWv027ntgiLERFJX14jv/8vgNvdfWu44V8rMxsPjAfo3bt3I5d0cFIPy5x22JEc9/ESCqsqwx26/xFNYSIi9ZDOlv4qoFdSv2c4VusyZpYHdADWAycDvzGzZcAPgP8yswmpH+Du97h7qbuXFhUV1ftLRGHSgC/UtEfNeS7CSkRE0pdO6M8E+ppZiZnlA6OBqSnLTAXGhe2LgRc88Fl3L3b3YuAO4CZ3vytDtUfqieM+z87c4A+lE9csgnnx2GUhIvFWZ+iHc/QTgKeBBcDj7j7PzG4wsxHhYvcSzOGXA1cB+xzWGTeftG7Ps32SrrN///3RFSMikqa05vTdfTowPWXsuqT2DmBUHe/xi4Oor1mbNPCLnLvo1aDz0EPw618HJ3CJiDRTOiO3AV4pPoE1bTsHnY8/hunTD/wCEZGIKfQboDonlykDhiUG/vzn6IoREUmDQr+BHhv4JaoJD0d96ikoL4+2IBGRA1DoN9DKjofz4lGlQccd/vjHaAsSETkAhX4GPDg46f65990H23SGrog0Twr9DHi55ESWHtot6HzyCTz8cLQFiYjsh0I/A9xyeOjEpK39u+7SjdNFpFlS6GfI5AHD2N6qIOjMnQuvvBJtQSIitVDoZ8jmwrb847j/lxi4447oihER2Q+FfgY9kLxD94knYNGi6IoREamFQv8AiidO2+eSygeyuKgYzj476LjDrbc2+meKiNSHQj/TJiZda+6hh2D16uhqERFJodDPtNNPh1NPDdqVlZrbF5FmRaGfaWbwk58k+n/6U3DsvohIM6DQbwznnQfHHhu0t2yBP/wh2npEREIK/caQkwNXX53o33YbbN4cXT0iIiGFfmMZOxZKSoL2+vVw553R1iMiQpp3zpL62XPI5bKf/xy+/vVg8Le/ZdC6vmxq3Y5lN58TXXEiktW0pd+Yxo6Ffv2C9ubNfGvmP6KtR0SynkK/MeXlwfXX13S/UTaVztt0JI+IREehn2H7nE07ahQMGADAIbt28N03JkVQlYhIQKHf2HJy4MYba7pfmz1Nt1QUkcgo9BvRnuvoFL+ey8we/QHIr66CH/844spEJFsp9JuCGTcO+2ai/8QT8OKL0dUjIllLod9E5nQ7minJ19u/8krYvTu6gkQkKyn0m9CtZ4xL3F3r3XeDm6iLiDQhhX4T+rB9F/508sWJgWuugXXroitIRLKOQr+J3TP0AujdO+isXw8/+lG0BYlIVlHoN7EdrQrh7rsTAw88AC+8EF1BIpJVFPoRKH7V4OKkaZ7LL4cdO6IrSESyhkI/Kr/7HXToELTffx9++cto6xGRrKDQj0q3bnDzzYn+zTfDzJnR1SMiWUGhH6Xx44N76kJwzP4ll1C4S9M8ItJ40gp9MxtuZovMrNzMJtbyfIGZPRY+/6aZFYfjXzSzWWY2N/z3zMyW38Ll5MCDD7I1v3XQX7yYa166P9qaRCTW6gx9M8sF7gbOAvoDY8ysf8pilwEb3b0PcDtwSzi+DjjP3QcA44CHMlV4XBT/eT7XDxtf0x83exqfWzKr5ro9IiKZlM6W/lCg3N2XuHsl8CgwMmWZkcADYXsyMMzMzN3fdvfV4fg8oLWZFWSi8DiZNOALPN33lJr+b6ffTtHWDRFWJCJxlU7o9wBWJvUrwrFal3H3KmAT0DllmYuA2e6+M/UDzGy8mZWZWdnatWvTrT0+zLhm+Pega1cAirZ9wl1Tf0Nuta7NIyKZ1SQ7cs3sOIIpn2/X9ry73+Pupe5eWlRU1BQlNTsb2nSAhx9mtwWr5OSV7/Hjlx+MuCoRiZt0Qn8V0Cup3zMcq3UZM8sDOgDrw35P4B/Ape7+QUMLjrVhw7jt9LE13cvfnAL//GeEBYlI3KQT+jOBvmZWYmb5wGhgasoyUwl21AJcDLzg7m5mHYFpwER3fy1TRcfZH04dxfNHDUkMfO1r8N57NV3t4BWRhqgz9MM5+gnA08AC4HF3n2dmN5jZiHCxe4HOZlYOXAXsOaxzAtAHuM7M3gkfh2X8W8SIWw5XnXMVKzsE8/ts2QLnngsffRRtYSISC3npLOTu04HpKWPXJbV3AKNqed0vAV1fIE17tuA3tW7Hty78KZMfvpq2lZ/C8uUwcqTutiUiDaYzcpuphYeVMGHE1cEJXABvvgnjxpGjI3pEpAEU+s3YS0cNgdtvTwxMmsSNz/4R3KMrSkRatLSmdyTz0t4Z+73vQXk5/P73AIx95ym2FBwCfg6YNWKFIhJH2tJv7szgjjvgkktqhi5/cwr8+tcRFiUiLZVCvyXIyQluoj4y6eoX116ra/CLSL0p9FuKVq3g0Ud59YhBibGf/SwIf83xi0iaFPotSWEh37zoZ7xcfGJi7Kab4Morobo6urpEpMVQ6DdzqTt8d7Qq5FsX/Sw4YWuPO++EMWNgx44Gna2rs31F4k+h3wLtzMuHKVOY1u+0xODjj8MXvkDHTzdHV5iINHsK/ZYqP5/vjbia+086LzH22mv8/aEfwfz5+yyuLXgRAYV+i1adk8v1w8Zz45nfrDlm/8iNq2HoUJg8OeLqRKQ5Uui3dGbcO+R8mDQJ2rQJxrZtg1Gj4OqrYdeuaOsTkWZFod+CHHCK5qKLYMYMlnXslhi79VY47bTgjN5030dEYk2hHycDBzJi3O1w9tmJsZkz4YQTGDXnGR3PLyIK/ZZqf1vrmwvbwr/+Bb/5TXBCF8C2bdz6798FPwbLl+/zPgfa8q/rOf3VINKyKPTjKCcHfvxjeOMN6NcvMf7UU3DccYyb9S/YrUs0i2QjhX4LUN8t6prlBw+GWbPg+9+nmvCKnNu2cf1zf2Z+9778x9hb9nqNiMSfQj/uDjkE7ryTiy/5DRx7bM1w/4+X8vjfJvL7f97CERtXR1igiDQlhX6WmN3jWHj7bbjxRj7NK6gZP2/hKzz//y/n1//+Hd03fxxhhSLSFHQTlWxSUAA//Slnru7JNS/dz4gFLwOQ59WMmfMMF857gUcGDWfolvW8dfelNS9LnvpZdvM5TV62iGSOtvSz0Jr2RXx/xNVcOPZWXjtiYM14we4qvj77SV7902Vw6aXBXwYiEisK/Zg70E7g2T2PZezomxgz+lfM6n5MzXh+dRU89BAMHswjj1zDsPI3074he1PtENaOZ5GDo+kdYcYRg7jokoF8fkkZ331jEkMrEhdsO3XFXE5dMZc1bTszZcAwKO8HffpEWK2INIRCXwJmvHTUEF46agiDVi/in7tmBtfzCY/n77Z1PRNmPA59H4czzoBLL6Xjp234pHX7iAsXkfpQ6MdY6hRIulMi73bvRzH96D7+S1w6exoXz32OLts3JRZ4+WV4+WXKLIfXjxgER66CCy6AoqK0atqzMzi5LSJNQ3P6sl+r2x/GzZ//Bqd89wHGX3Atz/YZCrm5Nc/neTVnLHsbvv1tOPxw+Mxn+N5rj0BZmW7fKNJMKfRjojGvg1OVm8czR5/Kty66Dioq4L//G045Ze+Fqqthxgx++OrDMGQI69p1grFjufbLV/DFy/5w0D8CzXWHra47JC2Vpnekfg4/HK66Cq66ilO/+z8MX/w6Zy98jSGrF+x1Fc8u2zfB3/7Gr/YMPPEzOO00xm/pAq8fCoMGRVK+SLZT6MtBW9O+iPtLR3J/6UiW/ehkeOYZptz0Fz63dPbe+wAANmyAf/2L/wI47X7IyeGZQ3vCqjOCawQNHgwnnAAdOkTxVUSyhkJfMqNLF/jqV/nhnA6YV9P/oyVMG1TNU3+aRGnF/H1/BKqrOXr9CvjrX4NH6MO2nTj85BOD6wQdcwynLt8Eq0+Ebt1qbgkpIgdPoS8Z55bDvMP7wJXncPlHR4M7R22ooLRiPkMr5nHR7jWwcGGtN3U5fOsGeP754AE8AvDotVBYCMXFwaOkhPELdsCk7VBSAr17Bz86OdpFJVIXhb7Uy56dl/U61NKMDzr34oPOvXhs0Jf5IdB6+A6OXbuUv5/SGmbPZv6TL3LU+pUU7K6q/T127Ah+KBYuBAimiV66P/F8Xh507co71W1Ye0gnvjjsBOjePfgLoaiIUZMWsaF1B56/6ULo1Any8nTIqGSltELfzIYDdwK5wF/c/eaU5wuAB4GTgPXAV9x9WfjcNcBlwG7g++7+dMaqlxbr0/zC4MqfVwShe/bEaeRW76bXJx/y0rmHc/Nt/6DP+pX0Wb+SE3ashY0bD/yGVVWwahUn7OmXv7nX05P2NO79TvBvx468RGt4sVfwV0LHjtC+fd2Pdu2gTRvyq3ZRmattJml5zOu4b6qZ5QKLgS8CFcBMYIy7z09a5rvAQHe/3MxGAxe4+1fMrD/BX+hDge7Ac8DR7r7fC7mUlpZ6WVnZQX8hHUYXP8tuPgc2bYJlyxh//SR6bvqYXps+5BvdgaVLYdWqun8UGsFuyyG3TWto0wZatw4eYfuVVdv57MBeibHCQsjPh/x87nptBVU5ufzg7OODsVatap5L7V/64Gx25bbikQlnBH/N5Obu+9jf+P4eOTnaPxJDZjbL3UvrWi6dTZWhQLm7Lwnf+FFgJDA/aZmRwC/C9mTgLjOzcPxRd98JLDWz8vD9ZqT7RUSA4KieQYN45uiKmqFvJE/N7NgBa9Zw4c+mcNi2DXTduoHrh3SCNWtg3TpmzS6n444tHMWnwQ9EBm4Sn+vVsG1b8EjxWYDl79T6ugl7Gq/V/RkP7mk8Wv/6DignZ98fA7O9H3t+HBpz7EDje6T+QKX7XH2WbU6fccUVcE7jTTumE/o9gJVJ/Qrg5P0t4+5VZrYJ6ByOv5Hy2h4HXa3I/hQWQkkJs3sm7g52/U2J/3EuSt4XsXs3fPIJZ06cTMdPt3Dojs3ce/7R/PTB12lXuZ22O7dzxUmHwebN+zzWrl5Lwa5KCqsqg6uRtlTV1cFj166oK5FUI0Y06ts3i0lJMxsPjA+7W81sUQPerguwruFVtSix/s52S63DXeyWA3/n2l63n/fivil79ye8UftyEYv1et6P7PvO3/lOF77znYP5zkeks1A6ob8K6JXU7xmO1bZMhZnlAR0Iduim81rc/R7gnnQKrouZlaUzrxUn+s7ZQd85OzT2d07nwOaZQF8zKzGzfGA0MDVlmanAuLB9MfCCB3uIpwKjzazAzEqAvsBbmSldRETqq84t/XCOfgLwNMEhm/e5+zwzuwEoc/epwL3AQ+GO2g0EPwyEyz1OsNO3CrjiQEfuiIhI40prTt/dpwPTU8auS2rvAEbt57W/gsR1t5pARqaJWhh95+yg75wdGvU713mcvoiIxIcuViIikkViE/pmNtzMFplZuZlNjLqexmBmvczsRTObb2bzzOw/w/FOZvasmb0f/nto1LVmmpnlmtnbZvZk2C8xszfD9f1YeJBBrJhZRzObbGYLzWyBmZ0a93VtZleG/22/Z2aPmFlh3Na1md1nZh+b2XtJY7WuVwv8Lvzuc8xscEM/PxahH14q4m7gLKA/MCa8BETcVAE/dPf+wCnAFeH3nAg87+59gefDftz8J7AgqX8LcLu79wE2ElzfKW7uBJ5y92OAQQTfP7br2sx6AN8HSt39eIIDR0YTv3X9P8DwlLH9rdezCI567EtwLtMfG/rhsQh9ki4V4e6VBCetj4y4poxz9zXuPjtsbyEIgR4E3/WBcLEHgPOjqbBxmFlP4BzgL2HfgDMJLvkB8fzOHYAzCI6Mw90r3f0TYr6uCQ4uaR2e79MGWEPM1rW7v0xwlGOy/a3XkcCDHngD6Ghm3Rry+XEJ/douFRHryz2YWTFwIvAm0NXd14RPfQh0jaisxnIHcDWw50a7nYFP3H3PdRDiuL5LgLXA/eG01l/M7BBivK7dfRXwW2AFQdhvAmYR/3UN+1+vGc+2uIR+VjGztsAU4Afuvjn5ufCkuNgckmVm5wIfu/usqGtpYnnAYOCP7n4isI2UqZwYrutDCbZsSwiuynsI+06DxF5jr9e4hH5al3uIAzNrRRD4D7v738Phj/b8yRf++3FU9TWC04ARZraMYNruTIK57o7hFADEc31XABXuvufGAJMJfgTivK6/ACx197Xuvgv4O8H6j/u6hv2v14xnW1xCP51LRbR44Vz2vcACd78t6anky2CMA/7Z1LU1Fne/xt17unsxwXp9wd3HAi8SXPIDYvadAdz9Q2ClmfULh4YRnNke23VNMK1zipm1Cf9b3/OdY72uQ/tbr1OBS8OjeE4BNiVNAx0cd4/FAzib4GYvHwDXRl1PI33H0wn+7JsDvBM+ziaY434eeJ/gRjWdoq61kb7/54Enw/aRBNdxKie4MVZB1PU1wvc9ASgL1/cTwKFxX9fA9cBC4D3gIaAgbusUQ0sFAAAAVklEQVSa4MZSa4BdBH/RXba/9QoYwZGJHwBzCY5satDn64xcEZEsEpfpHRERSYNCX0Qkiyj0RUSyiEJfRCSLKPRFRLKIQl9EJIso9EVEsohCX0Qki/wf9pVvWDlrSDUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 10000\n",
    "\n",
    "xmax = 60\n",
    "ymax = xlambda\n",
    "\n",
    "rx = np.random.random(N)*xmax\n",
    "ry = np.random.random(N)*ymax\n",
    "\n",
    "values = []\n",
    "\n",
    "Nin = 0\n",
    "for i in range(N):\n",
    "    if(ry[i] <= xlambda*np.exp(-xlambda*rx[i])):\n",
    "        # Accept\n",
    "        values.append(rx[i])\n",
    "        Nin += 1\n",
    "        \n",
    "x = np.asarray(values)\n",
    "\n",
    "print(\"Acceptance Ratio: \",Nin/float(N))\n",
    "\n",
    "binwidth=xlambda*5\n",
    "pyplot.hist(x,bins=np.arange(0.,100., binwidth),density=True);\n",
    "pyplot.plot(np.arange(0.,100.,binwidth),xlambda*np.exp(-xlambda*np.arange(0.,100.,binwidth)),ls='-',c='red',lw=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Challenge 9.1:\n",
    "\n",
    "Improve the acceptance ratio by using a linear function $f(x)=1-\\alpha x$, with a ppropriate choice of $\\alpha$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Random walk methods: the Metropolis algorithm\n",
    "---------------------------------------------\n",
    "\n",
    "Suppose that we want to generate random variables according to an\n",
    "arbitrary probability density $P(x)$. The Metropolis algorithm produces\n",
    "a “random walk” of points $\\{x_i\\}$ whose asymptotic probability\n",
    "approaches $P(x)$ after a large number of steps. The random walk is\n",
    "defined by a “transition probability” $w(x_i \\rightarrow x_j)$ for one\n",
    "value $x_i$ to another $x_j$ in order that the distribution of points\n",
    "$x_0$, $x_1$, $x_2$, ... converges to $P(x)$. In can be shown that it is\n",
    "sufficient (but not necessary) to satisfy the “detailed balance”\n",
    "condition $$p(x_i)w(x_i \\rightarrow x_j) = p(x_j)w(x_j \\rightarrow x_i).\n",
    "$$ This relation dos not specify $w(x_i \\rightarrow x_j)$\n",
    "uniquely. A simple choice is\n",
    "$$w(x_i \\rightarrow x_j)=\\min{\\left[ 1,\\frac{P(x_j)}{P(x_i)} \\right] }.$$\n",
    "This choice can be described by the following steps. Suppose that the\n",
    "“random walker” is a position $x_n$. To generate $x_{n+1}$ we\n",
    "\n",
    "1.  choose a trial position $x_t=x_n+\\delta _n$ , where the $\\delta _n$\n",
    "    is a random number in the interval $[-\\delta ,\\delta]$.\n",
    "\n",
    "2.  Calculate $w=P(x_t)/P(x_n)$.\n",
    "\n",
    "3.  If $w \\geq 1$ we accept the change and let $x_{n+1}=x_t$.\n",
    "\n",
    "4.  If $w \\leq 1$, generate a random number $r$.\n",
    "\n",
    "5.  If $r \\leq w$, accept the change and let $x_{n+1} = x_t$.\n",
    "\n",
    "6.  If the trial change is not accepted, the let $x_{n+1}=x_n$.\n",
    "\n",
    "It is necessary to sample a number of points of the random walk before\n",
    "the asymptotic probability $P(x)$ is attained. How do we choose the\n",
    "“step size” $\\delta$? If $\\delta$ is too large, only a small fraction of\n",
    "changes will be accepted and the sampling will be inefficient. If\n",
    "$\\delta$ is too small, a large number will be accepted, but it would\n",
    "take too long to sample $P(x)$ over the whole interval of interest.\n",
    "Ideally, we want at least 1/3-1/2 of the trial steps to be accepted. We\n",
    "also want to choose $x_0$ such that the distribution $\\{x_i\\}$ converges\n",
    "to $P(x)$ as quickly as possible. An obvious choice is to begin the\n",
    "random walk at the point where $P(x)$ is maximum.\n",
    "\n",
    "### Exercise 9.1: The Gaussian distribution\n",
    "\n",
    "1.  Use the Metropolis algorithm to generate a Gaussian distribution\n",
    "    $P(x)=A \\exp{(-x^2/2\\sigma ^2)}$. Is the numerical value of the\n",
    "    normalization constant $A$ relevant? Determine the qualitative\n",
    "    dependence of the acceptance ratio and the equilibrium time on the\n",
    "    maximum step size $\\delta$. One possible criterion for equilibrium\n",
    "    is that $\\langle x^2\n",
    "    \\rangle \\approx \\sigma ^2$. For $\\sigma = 1$, what is a reasonable\n",
    "    choice of $\\delta$? (choose $x_0 = 0$.)\n",
    "\n",
    "2.  Plot the asymptotic probability distribution generated by the\n",
    "    Metropolis algorithm.\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VdW58PHfQ0ICgoyJiARIlIAMMmgEcQSnQlWoFSvaWm1pqR/rfdvXTtheqdf29uqtV9perYpai9aKviiKgiCDzDKEeSoawhQEDVMYQwh53j/2zjk7x4ScwDlnn+H5fj7nw17rrL3PswmcJ3uvtdcSVcUYY4xp5HcAxhhj4oMlBGOMMYAlBGOMMS5LCMYYYwBLCMYYY1yWEIwxxgCWEIwxxrgsIRhjjAEsIRhjjHGl+x1AQ2RlZWlubq7fYRhjTEJZsWLFXlXNrq9dQiWE3NxcCgsL/Q7DGGMSiohsD6ed3TIyxhgDWEIwxhjjsoRgjDEGsIRgjDHGZQnBGGMMYAnBGGOMyxKCMcYYwBKCMcYYlyUEY4wxQII9qWxMvMsdMzWwve2JW3yMxJiGsysEY4wxQJgJQUSGiMhmESkSkTG1vJ8pIm+67y8VkdyQ9zuJyBER+Xm4xzQmXuWOmVrjZUyyqDchiEga8CwwFOgB3C0iPUKajQIOqGoXYBzwZMj7TwMfNvCYxhhjYiicPoT+QJGqFgOIyERgOLDR02Y48Ji7PQl4RkREVVVEvgFsBY428JjGJISLf/Y2g4pX0Hf3p3ynZTtm5A+ktHlrv8MypsHCSQgdgJ2ecgkwoK42qlopImVAWxEpB34F3AT8vLb2pzkmACIyGhgN0KlTpzDCNSY2Lv5yK/+2eCKDiws55+SJQP3jHz3H8pweUHAcRozwMUJjGibancqPAeNU9ciZHkBVx6tqgaoWZGfXu76DMTExYMc63v7HL7hl86IayQCgEcqAkg1w553wZOjdU2PiVzhXCLuAjp5yjltXW5sSEUkHWgL7cH7rHyEi/w20Aqrcq4YVYRzTmLjh7Ty+eusqXnzn9zStDCaCT9t2Yn5eP3p+WUz/nRtI0yrnjTFjoLwcxo4FkViHbUyDhJMQlgP5IpKH86U9ErgnpM0U4D7gE2AEMEdVFbimuoGIPAYcUdVn3KRR3zGNiTuDthTywuT/JPPUSQC+aN6GH3zzUda1zw+0yTp6gD+//0eu2r7WqXjsMf46fT3/fd39gD2fYOJXvbeMVLUSeAiYAWwC3lLVDSLyuIgMc5u9jNNnUAQ8DJx2GGldxzzz0zAm+tofKuWZKU8GksGuc7P51j1P1EgGAHubteb7d/yWuXmXBeoeXDKJr/9rYUzjNaahxPlFPjEUFBSorals/JD7qw948Z3fcVPRMgBKWpzHyHv+i5KW7ercJ6PyJH999w/cuGU5AF82a82NP3iOtX+6KyYxG1NNRFaoakF97exJZWPCMOTTxYFkAPDT23522mQAUJHemIdv/RlfNG8DwHlHD/CreX+PZpjGnBVLCMbU5+BB/mPWC4HiP/sMoTCnZ1i7HmrSnN/e+KNA+durp8NCu3Vk4pMlBGPqUD01xT8G30O7I/sB57bPE4Pub9Bxpne9kpldPI/ZjB4NJ07UvYMxPrGEYMxp5O3fxT2rpwfKj934Iw41ad6wg4gw9qYHOJLR1Clv2sS/D3vY5kEycccSgjGn8cNl79AIZ+DF/Nx+TOt21RkdZ3eLbP585d01jptWdSoiMRoTKZYQjKlD9pH93LF+dqD87MBvndXDZa/3G8pB9+qi88E9DN286KxjNCaSLCEYU4dRhe+ReaoSgJUXdGNpx15ndbxjGU2ZcOmtgfIDS9+GBBr2bZKfJQRjalNWxj2rAjO28/yAERGZemLCZbdxPD0TgF5fbIFZs876mMZEiiUEY2rz/PO0qDgGQFGbHGbm1zoZb4PtP6clb/a+KVhhk9+ZOGIJwRiP3DFT6frzd/nyd8Ev6hcGfBOVyP1Xean/7VRWH2/2bLCn702csIRgTIgbipZy3tEDAOxp3ob3egyO6PFLWrbjg+7XBCvGj4/o8Y05U5YQjAkxwjOyaGKfr1GR3jjin/FaP8+Mp2++CceORfwzjGkoSwjGeGQf2c+g4hWB8qReN0Tlc1Z06E5x6wucwqFD8O67UfkcYxrCEoIxHrdv+DiwuM0nnS6hpNX50fkgESZdcmOw/Mor0fkcYxrAEoIx1VS5c11wGOikXjeepvHZm9xzcHAo6+zZsGNHVD/PmPpYQjCm2vLl5O/bCcCRjKZnPE1FuHa3yGZ+575OQRVefTWqn2dMfcJKCCIyREQ2i0iRiHxlNTQRyRSRN933l4pIrlvfX0RWu681InK7Z59tIrLOfc/G3Rn/eW7bTOt2FcczmkT9I2vcNvr73+3JZeOrehOCiKQBzwJDgR7A3SLSI6TZKOCAqnYBxgHVg7jXAwWq2hcYArzgrqdcbbCq9g1nJR9joqq8HN54I1Cs8UUdRTPyr+BQZjOnsGWLrZVgfBXOFUJ/oEhVi1W1ApgIDA9pMxyY4G5PAm4QEVHVY+76yQBNAPv1x8SnDz+EsjIAtrc6n2VhLoBztk40zuR97zMJnqRkTKyFkxA6ADs95RK3rtY2bgIoA9oCiMgAEdkArAMe8CQIBT4SkRUiMvrMT8GYCHj77cDmlO7XRWTeonBN6X5dsDB5MlRVxeyzjfGKeqeyqi5V1Z7A5cAjIlJ9Y/ZqVb0U51bUj0Xk2tr2F5HRIlIoIoWlpaXRDtekohMn4P33A8Xp3a6M6ccvz+lB6TmtnMKePbB4cUw/35hq4SSEXUBHTznHrau1jdtH0BLY522gqpuAI0Avt7zL/fNLYDLOramvUNXxqlqgqgXZ2dlhhGtMA82e7TwcBuxo2Y4N510Y04+vapTGzPwrghWeqxVjYimchLAcyBeRPBHJAEYCU0LaTAHuc7dHAHNUVd190gFEpDNwMbBNRJqJyLlufTPgZpwOaGNiz/MF/GG3q2J6uyj4uZ6rknfesdFGxhf1JgT3nv9DwAxgE/CWqm4QkcdFZJjb7GWgrYgUAQ8D1UNTrwbWiMhqnKuAB1V1L9AOWCgia4BlwFRVDS5ca0ysVFbCe+8FitO7xvZ2UbVPOvWG1q2dwo4dNgOq8UV6/U1AVacB00Lqxnq2y4E7a9nvNeC1WuqLgT4NDdaYiJs3D/a5dzc7dGD1BV19CaMyLR2GDYMJ7mC9t9+Gyy/3JRaTuuxJZZPaPLeLXjn/0oiue9BQo450ChbetuU1TexZQjCpq6rKGebpmtF1oI/BwMLcfhzJaOoUiopgvXWrmdiyhGBS15IlzjBPYO85LWP2MFpdTqRn8PGFnof2bbSRiTFLCCZ1ffBBYHNWlwFUNUrzMRhHjU7tqVP9C8SkJEsIJnV5vnDnXBQfHbgL8voF11suLAxcwRgTC5YQTGoqKYG1awGoaJTOwty+PgfkONSkOYU5nrkjp9tobBM7lhBMapoWHEW9tGMvjlV35sYB79XK1Cde8jESk2osIZiU9NFTwbUPPo6T20XVvB3L12xdBSdP+hiNSSWWEEzqOXGCq7atDhQ/vii+luP4LKsTJS3OA6BFxTFYtMjniEyqsIRgUs/8+TQ7WQ7A1tbt2domdDZ3n4nUTFLTptXd1pgIsoRgUo9ndNHcC+Pr6qBajVFPNvzUxIglBJN6PL9xx8tw01CfdLqE8vQMp7BxI2zb5ms8JjVYQjCp5bPPnBdwrHEmyzr28jmg2pU3bsInnS4JVthtIxMDlhBMavnww8Dmos59OFH9W3gcqnH14onbmGixhGBSy8yZgc15eZf5GEj95uddGizMnWvDT03UWUIwqePkSeeL1bUgr59/sYRhe6v27GzZzikcOeJMxmdMFIWVEERkiIhsFpEiERlTy/uZIvKm+/5SEcl16/uLyGr3tUZEbg/3mMZE3JIlzhcrsLNlO7a3au9zQPUQYUGuJ2l99JF/sZiUUG9CEJE04FlgKNADuFtEeoQ0GwUcUNUuwDjgSbd+PVCgqn2BIcALIpIe5jGNiSzPF+qC3L6+rJ3cUAs8cyytevn/+RiJSQXhXCH0B4pUtVhVK4CJwPCQNsMBd+0/JgE3iIio6jF3TWaAJkD1ElDhHNOYyPL0HyzMje/bRdUWd+5DFU7i6r3nMzhwwOeITDILJyF0AHZ6yiVuXa1t3ARQBrQFEJEBIrIBWAc84L4fzjGNiZwDB2D5cmdbhEWdE2NJ77Km57K2fRcA0rQK5szxOSKTzKLeqayqS1W1J3A58IiINGnI/iIyWkQKRaSwtLQ0OkGa5Pfxx86SmQAFBZQ1PdffeBqgxtWM5yrHmEgLJyHsAjp6yjluXa1tRCQdaAns8zZQ1U3AEaBXmMes3m+8qhaoakF2dnYY4RpTC2+H7E03+RfHGVhgCcHESDgJYTmQLyJ5IpIBjASmhLSZAtznbo8A5qiquvukA4hIZ+BiYFuYxzQmcrxfpDff7F8cZ2Blh4s52ti9sC4uhi1b/A3IJK16E4J7z/8hYAawCXhLVTeIyOMiMsxt9jLQVkSKgIeB6mGkVwNrRGQ1MBl4UFX31nXMSJ6YMQHFxc4LoFkzGDjQ33ga6GRaY5Z6p9iwqwQTJenhNFLVacC0kLqxnu1y4M5a9nsNeC3cYxoTFZ4v0DntuvP9sYn3hbowtx/XFxc6hZkz4YEH/A3IJCV7UtkkP8/InEVxsnZyQ3mfR6jRQW5MBFlCMMmtqsr5AnUlynDTUJ9ldaL0nFZO4cABWLPG34BMUrKEYJLb+vXgDlfe17QFm7M7+xzQGRLhk869g+XZs/2LxSQtSwgmuXluF33SuQ8qiftPvsbVjT2gZqIgcf93GBOGmc++Edhe7P0NOwEt9iaE+fNtOmwTcZYQTPKqrGTAjvWB4uJOiZ0QdrZsB53dW15Hjwan4jAmQiwhmOS1ciUtKo4B8Pm5WWxrfYHPAZ0lEbj++mDZbhuZCLOEYJKXp+P1k869E2K66/r8tLRNsGAdyybCLCGY5OV9/iBBh5uGqnHba/FiOH7cv2BM0rGEYJLTiROwcGGg+EmC9x9U+/LcthS1yXEKFRVOUjAmQiwhmOS0ZAmUlwNQ3PoCdrdInplyF9vwUxMllhBMcqrx/EFyXB1UW2wPqJkosYRgkpM3ISTJ7aJqSzpdEiwUFsLhw/4FY5KKJQSTfI4dg6VLA8UaX6BJ4GDTFmzKznUKp07V6Csx5mxYQjDJZ/HiwFO8n7btxN5mrX0OKPJqXPV4Ju8z5mxYQjDJx/MF+Unn5Lo6qFajX8QSgomQsBKCiAwRkc0iUiQiY2p5P1NE3nTfXyoiuW79TSKyQkTWuX9e79lnrnvM1e7rvEidlElx3oSQZP0H1ZZ27BV80G7lSigr8zcgkxTqTQgikgY8CwwFegB3i0iPkGajgAOq2gUYBzzp1u8FblPVS3DWXA5dPe3bqtrXfX15FudhjOPIkRpz/NRYejKJHGrSHPr1cwpVVbBggb8BmaQQzhVCf6BIVYtVtQKYCAwPaTMcmOBuTwJuEBFR1VWq+rlbvwFoKiKZkQjcmFotXAiVlc52794cOKelv/FE0fg0z9oOdtvIREA4CaEDsNNTLnHram2jqpVAGdA2pM0dwEpVPeGpe8W9XfSoSBJMNGP85/1iHDTItzBiocboKUsIJgJi0qksIj1xbiP9yFP9bfdW0jXu69469h0tIoUiUljqrnxlTJ28X4yDB/sXRwws79iTU9UL/qxeDfv3+xuQSXjhJIRdQEdPOcetq7WNiKQDLYF9bjkHmAx8V1W3VO+gqrvcPw8D/8S5NfUVqjpeVQtUtSA7O3mmHzBRcOgQrFjhbIvAddf5G0+UHc5sxrrzL3IKqs6iOcachXASwnIgX0TyRCQDGAlMCWkzBafTGGAEMEdVVURaAVOBMaq6qLqxiKSLSJa73Ri4FViPMWdjwQKngxVYf96F5D6Z/BO/LfGOopo717c4THKoNyG4fQIPATOATcBbqrpBRB4XkWFus5eBtiJSBDwMVA9NfQjoAowNGV6aCcwQkbXAapwrjBcjeWImBdUYbpqczx+EsgfUTCSlh9NIVacB00Lqxnq2y4E7a9nv98Dv6zjsZeGHaUwYPF+IyTZdRV2W5/SgUhqRrlWwdi3s3QtZWX6HZRKUPalsksOBA7BqFQCnpBHLkvT5g1DHMpqytn1+sGLePP+CMQnPEoJJDgsWOB2rwPp2F3E4s5nPAcWO97bRhD+84mMkJtFZQjDJIQX7D6p5E8IVO9b5GIlJdJYQTHKo0X+QnPMX1aUwpzsVjZzuwG57d8CXNguMOTOWEEzi27cP1qxxttPSWJ4TOtVWcitv3ITVF3QNVtjwU3OGLCGYxOd9IKuggKOZ5/gXi09s+KmJBEsIJvGl0HQVdVliCcFEgCUEk/gsIbDygm6cSHMfK9q8GT7//PQ7GFMLSwgmsZWWwnp31pP0dLjqKn/j8cmJxpmsuuDiYIU9j2DOgCUEk9g8X3yF7fLJ/d1c/2LxmfUjmLNlCcEkthRYLjNcts6yOVuWEExiS+EH0kKtbt+N4+nugoRFRbBz5+l3MCaEJQSTuHbvhk2bADiRls6KDt19DshfFemNKfT+HdhVgmkgSwgmcc2ZE9hc2aE7Jxrbct01bht5/n6MCYclBJO4PF94i1O8/6Da4s59goU5cwIT/hkTDksIJnF5EsKizn19DCR+rDu/C4czmjqFnTu57oGX/A3IJJSwEoKIDBGRzSJSJCJjank/U0TedN9fKiK5bv1NIrJCRNa5f17v2ecyt75IRP4iIhKpkzIpYOtW2LbN2W7WrOaaACnsVKM0lnrWgrhy+1ofozGJpt6EICJpwLPAUKAHcLeIhM4eNgo4oKpdgHHAk279XuA2Vb0EZ83l1zz7PAf8EMh3X0PO4jxMqvHeH7/2WirTwlr8LyV84rltdOX2NT5GYhJNOFcI/YEiVS1W1QpgIjA8pM1wYIK7PQm4QUREVVepavUz9BuApu7VRHughaouUVUFXgW+cdZnY1KHNyFcf33d7VLQYk/H8sAda60fwYQtnITQAfAOaC5x62pto6qVQBnQNqTNHcBKVT3hti+p55jG1E7VEsJp/Cs7l/1NWwCQdawMNmzwOSKTKGLSqSwiPXFuI/3oDPYdLSKFIlJYWloa+eBM4vnXv2DPHme7dWvo0+f07VOMSqOaD+nZ8FMTpnASwi6go6ec49bV2kZE0oGWwD63nANMBr6rqls87XPqOSYAqjpeVQtUtSA7OzuMcE3S837BDRoEaWm+hRKvPgkdfmpMGMJJCMuBfBHJE5EMYCQwJaTNFJxOY4ARwBxVVRFpBUwFxqjqourGqrobOCQiV7iji74LvHeW52JSxezZwW27XVSrGs9lzJ0LlZW+xWISR70Jwe0TeAiYAWwC3lLVDSLyuIgMc5u9DLQVkSLgYaB6aOpDQBdgrIisdl/nue89CLwEFAFbgA8jdVImiZ06VWNKhhtXp5E7ZqqPAcWn4jYd2NO8jVMoK4MVK/wNyCSEsMbqqeo0YFpI3VjPdjlwZy37/R74fR3HLAR61faeMXVasQIOHgRgT/M2FLXtWM8OKUqERbl9uWO9e7to1iwYMMDfmEzcsyeVTWKZOTOwuSi3L9jzjHVakNsvWPD8vRlTF0sIJrHMmhXYXGjTVZzWIm/H8uLFcPSof8GYhGAJwSSOo0dhUWBsAgtzLSGcTmnzNvwrq7NTOHkS5s/3NyAT9ywhmMSxYIHzxQbQsyel1Z2mpk6LPEnzxd+O9zESkwgsIZjE4b0PftNN/sWRQBZ4EsLV21b7GIlJBJYQTOLw9B98r6Slj4EkjmUde1HRyBlM2L10W/AJb2NqYQnBJIYvvoC1zlTOJ0OmeDZ1O5bRlFUdLg5WeB/qMyaEJQSTGDxfZCsvuJhj1YvAmHot9I428lxlGRPKEoJJDJ7+Axtd1DALQ59HsOmwTR0sIZj4p/rVB9JM2Na2z+dQZjOnsGsXbNrkb0AmbllCMPFv/Xrniwwoy2zGmvZdfQ4osZxqlFbzIbXp0/0LxsQ1Swgm/nm+wBbk9uNUI5vuuqHm5V0aLFhCMHWwhGDin+cLbN6Fl/kYSOKq8fc2b55NY2FqZQnBxLcjR5wnlF01ftM1YdvdIpvNWZ2cQkUF99//R38DMnHJEoKJbx9/HJyuondvvjw3dKluE655ecGrhOu2rvQxEhOvLCGY+PahZ92kIUP8iyMJeG8bDSou9DESE68sIZj4pVozIQwd6l8sSWB5Tk+ONc4EIO/Abigq8jkiE2/CSggiMkRENotIkYiMqeX9TBF5031/qYjkuvVtReRjETkiIs+E7DPXPWbo0prGOD77DLZtc7abN4crr/Q1nERXkd645lrLM2b4F4yJS/UmBBFJA54FhgI9gLtFpEdIs1HAAVXtAowDnnTry4FHgZ/Xcfhvq2pf9/XlmZyASWKe0UUfte9F7lhb9ets1RhtZMNPTYhwrhD6A0WqWqyqFcBEYHhIm+HABHd7EnCDiIiqHlXVhTiJwZiG8Q43tdFFEeHtWGbOHDhxwr9gTNwJJyF0AHZ6yiVuXa1tVLUSKAPCGQ7yinu76FGR2hfHFZHRIlIoIoWlpaVhHNIkhaNHnRFGLnv+IDJ2tG5PcesLnMKxY84zCca4/OxU/raqXgJc477ura2Rqo5X1QJVLcjOzo5pgMZHs2ZBuXNhuTmrEyUt2/kcUPL4+KLLg4UpU/wLxMSdcBLCLqCjp5zj1tXaRkTSgZbAvtMdVFV3uX8eBv6Jc2vKGMf77wc2Z3exfxqRNMv79/n++zb7qQkIJyEsB/JFJE9EMoCRQOivFVOA+9ztEcAc1br/lYlIuohkuduNgVuB9Q0N3iSpqir44INAcVaXAT4Gk3yW5/QMzn66YwdDRj3rb0AmbtSbENw+gYeAGcAm4C1V3SAij4vIMLfZy0BbESkCHgYCQ1NFZBvwNHC/iJS4I5QygRkishZYjXOF8WLkTssktOXLnRXSALKzWW2zm0ZUZVo6cz19MjcWLfUxGhNP0sNppKrTgGkhdWM92+XAnXXsm1vHYa2X0NTOc7uIW26hymY3jbhZXfozbNN8AG4sWuZzNCZe2JPKJv54Ozpvu82/OJLY3AsLqBTnv3/f3Z/C7t0+R2TigSUEE1+2b4d165ztjAy4+WZ/40lSh5o0Z3nHnsGKqVP9C8bEDUsIJr54bxddf70zZYWJilkXhYw2MinPEoKJL3a7KGZqDOedOROOH/cvGBMXLCGY+HHwIMydGygOXH8OuWPsVka0bGvTgaI2OU7h+HH46CN/AzK+s4Rg4sf77wcWw1l7fhd2t7An06NtRteBwcLbb/sXiIkLlhBM/Jg0KbA5rdvVPgaSOqZ1uypYmDLFJrtLcZYQTHw4dKjG/PwfdrO1D2JhQ7uL2FE9T1RZmTOHlElZlhBMfPjgg+Bvp336sL16Rk4TXSI1rxI8V2km9VhCMHFh+n8EF9R7qsUlPkaSej70JoR334WKCv+CMb6yhGD8d+QIg4pXBIo1vqBM1K1p35Vd57od+AcP1liHwqQWSwjGf9Om0aTS+a10c1YntrTtWM8OJqJEavbZ2G2jlGUJwfjP8wVkVwf+qDGqa/JkqKz0LxjjG0sIxl9Hj9aYR2eaJQRfrOrQjT3N2ziFffvstlGKsoRg/PXuu87avkBRmxw+zersc0CpSaVRzauz11/3LxjjG0sIxl+vvRbYnNxzMIj4GExqm9xzcLDw9tuBRG1SR1gJQUSGiMhmESkSkTG1vJ8pIm+67y8VkVy3vq2IfCwiR0TkmZB9LhORde4+fxGxb4KUs2ePM6ma672eg/yLxbD2/Hy2VM9tdOQIvPeevwGZmKs3IYhIGvAsMBToAdztLoPpNQo4oKpdgHHAk259OfAo8PNaDv0c8EMg330NOZMTMAnsjTec9ZOBpR17UVL9xKzxhwiTvUnZc/VmUkM4Vwj9gSJVLVbVCmAiMDykzXBggrs9CbhBRERVj6rqQpzEECAi7YEWqrpEVRV4FfjG2ZyISUD/+Edgc3KPQf7FYQLe9d42+uij4NrWJiWEkxA6ADs95RK3rtY2qloJlAFt6zlmST3HBEBERotIoYgUlpaWhhGuSQgbN8LKlc52ZibTLrbJ7OJBSct2cM01TuHUKZg40d+ATEzFfaeyqo5X1QJVLcjOtumQk4bn6oDbbuNQE1sZLW585zvBbe/PySS9cBLCLsD76GiOW1drGxFJB1oC++o5Zk49xzTJqqqqxrDGH1R19zEY8xV33umsZw1QWAibNvkbj4mZcBLCciBfRPJEJAMYCUwJaTMFuM/dHgHMcfsGaqWqu4FDInKFO7rou4ANaUgVs2bBjh0AHGhyLvMuvNTngIxX7pOL+TC3IFjxt7/5F4yJqXoTgtsn8BAwA9gEvKWqG0TkcREZ5jZ7GWgrIkXAw0BgaKqIbAOeBu4XkRLPCKUHgZeAImAL8GFkTsnEvRdeCGy+0+t6TqY19jEYU5s3e98ULLzyCpSX193YJI30cBqp6jRgWkjdWM92OXBnHfvm1lFfCPQKN1CTJD7/vMb49n/2sdHG8Wh+3qWUtMgm51CpM5XFO+/APff4HZaJsrjvVDZJ5m9/c0avAEs69mJLls1sGo+qGqXxhjdZP/+8f8GYmLGEYGLn1CkYPz5Q/GffoT4GY+rzVu+bONkozSksWAAbNvgbkIk6SwgmdqZPh53uIy1ZWUzvausmx7PS5m2Y2WVAsMKTzE1ysoRgYsd72+F736Mi3TqT412Nq7gJE2zCuyRnCcHExrZtMC04LuG6snz/YjFhW5Tbh22t2juFsjJn/imTtCwhmNj4058CE9nNz+3H9tYX+ByQCYdKI173XiU8/XTg52iSjyUEE33798OLLwaKL11u8xgmkjf73MyRjKZOYePGGld6JrlYQjDR99xzgXvPm7JzmZ9nTyYnkkNNmvNGn68FK/74R/+CMVFlCcFEV3k5/O//Borj+3/TVkVLQH8rGB4cgjp/Pixb5m9AJiosIZjoeu214Jz6OTm83/1af+MxZ2R3i+yaPzu7SkhKlhBM9FQuHs/lAAANn0lEQVRVwVNPBcs//SmVaWHNlmLi0Iv9bw8W3nkHior8C8ZEhSUEEz2TJsGnnzrbLVvCD3/obzzmrGw670LmVff/VFXBf/2XvwGZiLOEYKKjshIefTRQfLb7zeT+YYGPAZlIeO6KEcHC3/8Omzf7FouJPEsIJjomTAhcHRzKbMYLA+7wOSATCUs69WZB575OoaoKxo49/Q4moVhCMJFXXg6PPRYoPj/gDlsiM4k8de29wcJbb8GqVf4FYyLKEoKJvOefh5ISAErPacUrlw2rZweTSNZc0I3pXQcGK/793/0LxkRUWAlBRIaIyGYRKRKRMbW8nykib7rvLxWRXM97j7j1m0Xka576bSKyTkRWi0hhJE7GxIHDh+E//zNQfObKuzie0cTHgEw0PHXNvVThPk8ybZozPbZJePUmBBFJA54FhgI9gLs9y2BWGwUcUNUuwDjgSXffHjhrMPcEhgB/dY9XbbCq9lXVAkxyePxx2LsXgJIW59VcZMUkjaKsTkzuNThY8ZOfBBY+MokrnCuE/kCRqharagUwERge0mY4MMHdngTcICLi1k9U1ROquhVn/eT+kQndxJ3162HcuEDxqWvvtSmuk9jTV3+H8vQMp7BqlTNFiUlo4SSEDsBOT7nErau1japWAmVA23r2VeAjEVkhIqMbHrqJK6osG3pXjeUx3+0xyN+YTFTtanke/zvwrmDFb34De/b4F5A5a352Kl+tqpfi3Ir6sYjUOqeBiIwWkUIRKSwtLY1thCZ8r75K/5KNAJxslMa/3/ygzVmUAl7s/03o1s0pHDoEP/+5vwGZsxJOQtgFeFdCz3Hram0jIulAS2Df6fZV1eo/vwQmU8etJFUdr6oFqlqQnZ0dRrgm5vbvh1/8IlB86fLbKcrq5GNAJlYq0htzT1/PMNTXX4c5c/wLyJyVcBLCciBfRPJEJAOnk3hKSJspwH3u9ghgjqqqWz/SHYWUB+QDy0SkmYicCyAizYCbgfVnfzom5lThRz8C9+pt17nZ/OXKkT4HZWJpcW5f3ut+XbDie9+Dgwf9C8icsXoTgtsn8BAwA9gEvKWqG0TkcRGpHmD+MtBWRIqAh4Ex7r4bgLeAjcB04MeqegpoBywUkTXAMmCqqk6P7KmZmHjlFWfOItfYmx+wYaYp6PfXj+Jg9cOHO3bAAw84vyyYhCKaQD+0goICLSy0RxbixqefQr9+gcVvXu87hN987SGfgzJ++drmxbzw7h+CFX//O9x3X53tTeyIyIpwhvfbk8rmzFRUwD33BJIBF1/M767/gb8xGV/N6HYl/MDzb+DHP7YpshOMJQTTcKrw0EOwYoVTzsiAN96gvLHdKkp5f/oTdO3qbB89Crff7ow+MgnBEoJpuKeeghdfDBR/d/V3yZ0YOvDMpKLc383lloEPUtHIXQhp/Xq46y5nOnQT9ywhmIZ55x341a8Cxbd7DublgtAH100q23B+F8YM/bdgxfTpztQWCdRfmaosIZjwLV4M3/lO4D/20o69eGTI/7EH0MxXvNPrBv7ifYr5r3+tuZyqiUuWEEx45s2Dm2+G48cB2Nq6PT+6/dc2V5Gp09PXfIcp3T0TEPzyl/Dkk/4FZOplCcHUb9YsGDrU6SQE9p7Tku+PeIyDTVv4HJiJayL84us/ZWnHXsG6MWOcGXHt9lFcsoRgTm/SJLj11sCVAe3bc9fdT7C1Tej8hsZ81Yn0DL434rcs7tQ7WPnb38LPfmYdzXHIEoKp3alT8MgjcOedcOIEAJ+fm8WgWx9jS1bHenY2JuhYRlO+P2Is83P7BSvHjXOuOvft8y8w8xWWEMxX7d0Lt9wCTzwRqNrWqj3fuucJttmVgTkD5Y2b8MM7HmVG/hXBylmzoKAAVq70LzBTgyUEE6QKb7wBPXrAjBnB+qFDGXbfOEpane9fbCbhnUjP4IHbf82fr7w7WLltG/Tv76ylUF7uW2zGYQnBOIqLYdgwZzoKz7oTfxl4Fxf2eoBD1ROXGXMWVBox7ppv84NvPsrhjKZO5alT8Ic/sKVDvk2d7TNLCKmupMSZmbJbN/jgg2B9Tg73j/gtT197L1WN0ure35gzMCt/ALfd96caI5Au2l8CN9wAN94IS5b4GF3qsoSQqgoLYdQo6NIFXnihxoiPV/vdQq8RTzP3ost9DNAku21tOjDy7j/wm5sfDF4tAMyeDQMHOs+9vPuujUaKIZv+OpWUlMDkyfDaa7B8+Vffv+YavtlpGCtzusc+NpPSzj+0l/+78HXuWD+bdK2q+WZOjjON9h13QN++9mT8GQh3+mtLCMns+HHn0nvuXKeTeOnSWputbp/P/1xzLwty+9l/NuOrvP27+OnCf3Lbpvk0opbvptxc57mYwYPh2mshKyvmMSaiiCYEERkC/BlIA15S1SdC3s8EXgUuw1lL+S5V3ea+9wgwCjgF/B9VnRHOMWtjCaEOqrB7tzP3/MaNsGqV81qzxlm3oBYn0tKZevE1/KPv11nZ4WJLBCau5Bzcwz1rpvOttTPJOlZWd8Nu3ZxFmvr1g0sugfx8J2mkp8cs1kQQsYQgImnAp8BNQAnOGst3q+pGT5sHgd6q+oCIjARuV9W7RKQH8AbQH7gAmAW4k6Wf/pi1SYmEoOoMvzt61HkdPgxlZc5r/35nBFBpKXzxhXMLqKQEtm8PTCtxOicbpfFJp97M6DqQad2u4sA5LWNwQsacuYzKkwwuXs6QzYu5oWgZLSqO1b9Tejp06uTcasrJgfbtITvbeWVlQcuWzqtFC2jWzHmdcw40St4u1XATQjhptD9QpKrF7oEnAsNx1kmuNhx4zN2eBDwjIuLWT1TVE8BWd83l/m67+o4ZGc89B1Onhte2ruTorVcNlqu3Q19VVcE/q1+nTjmvysrg6+RJ51VR4TwNXF4eeCo4Uora5LCk0yUs7diL+XmXUtb03Ige35hoqkhvzIyuVzKj65U0PnWSATvWM3DHWq7YsY7eez6jcdWpr+5UWekMoy4ubtiHNW4MTZpAZqaz6FNGhlOXnh58paUFX40aOVfWoX+GvqDmdnW5LnW9160b/M//NOycGiichNAB2OkplwAD6mqjqpUiUga0deuXhOxb/ahrfceMjPXrw08ICawssxlb21zA1tYd2HReHhvaXcTG8/LsKsAkjZNpjVmY14+Fec4UGE1OlnNx6XZ6frGFHl8Wc9G+EnIPfM75R/af4Qe4v6AdPhzBqCNo796of0Tc32gTkdHAaLd4REQ2+xnPGcoCovvTPHEUdn/mvDbOjepHhSn65xyfUvG8fTvnzcB7fnywI7bnvXTp2fT1dQ6nUTgJYRfgnc0sx62rrU2JiKQDLXE6l0+3b33HBEBVxwPjw4gzbolIYTj375JJKp4zpOZ5p+I5Q3Kedzi9KMuBfBHJE5EMYCQwJaTNFOA+d3sEMEed3uopwEgRyRSRPCAfWBbmMY0xxsRQvVcIbp/AQ8AMnCGif1PVDSLyOFCoqlOAl4HX3E7j/Thf8Ljt3sLpLK4EfqyqpwBqO2bkT88YY0y4EurBtEQlIqPdW18pIxXPGVLzvFPxnCE5z9sSgjHGGMAmtzPGGOOyhBBlIvIzEVERyXLLIiJ/EZEiEVkrIpf6HWMkicgfReRf7rlNFpFWnvcecc97s4h8zc84I01EhrjnVSQiY/yOJ1pEpKOIfCwiG0Vkg4j8xK1vIyIzReQz98/WfscaaSKSJiKrROQDt5wnIkvdn/mb7gCZhGYJIYpEpCNwM7DDUz0UZ7RVPs7zFc/5EFo0zQR6qWpvnOlJHgFwpzEZCfQEhgB/dadFSXjueTyL87PtAdztnm8yqgR+pqo9gCuAH7vnOgaYrar5wGy3nGx+AmzylJ8ExqlqF+AAzpxtCc0SQnSNA34JNaZtHA68qo4lQCsRae9LdFGgqh+pavUE9ktwnjEBzzQmqroV8E5jkugC07uoagVQPRVL0lHV3aq60t0+jPMF2QHnfCe4zSYA3/AnwugQkRzgFuAltyzA9ThT9UCSnLMlhCgRkeHALlVdE/JWbVOBJOvK9d8HPnS3k/m8k/nc6iQiuUA/YCnQTlV3u2/tAdr5FFa0/Annl7vqxRraAgc9v/wkxc887qeuiGciMguobeX53wC/xrldlHROd96q+p7b5jc4txdej2VsJjZEpDnwNvBTVT0knikVVFVFJGmGL4rIrcCXqrpCRAb5HU80WUI4C6p6Y231InIJkAescf+j5AArRaQ/4U0FEtfqOu9qInI/cCtwgwbHNSf8eZ9GMp/bV4hIY5xk8LqqvuNWfyEi7VV1t3sL9Ev/Ioy4q4BhIvJ1oAnQAmctl1Yiku5eJSTFz9xuGUWBqq5T1fNUNVdVc3EuJy9V1T04U3R81x1tdAVQ5rnUTnjuwke/BIapqnfy+rqmMUkGKTMVi3vv/GVgk6o+7XnLO33Nffg651xkqeojqprj/l8eiTM1z7eBj3Gm6oEkOWe7Qoi9acDXcTpVjwHf8zeciHsGyARmuldHS1T1gdNNY5Lo6prexeewouUq4F5gnYisdut+DTwBvCUio4DtwLd8ii+WfgVMFJHfA6twEmVCsyeVjTHGAHbLyBhjjMsSgjHGGMASgjHGGJclBGOMMYAlBGOMMS5LCMYYYwBLCMYYY1yWEIwxxgDw/wHUp72vlvJtiwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 100000\n",
    "\n",
    "x = np.zeros(N)\n",
    "delta = 2.\n",
    "sigma = 10.\n",
    "sigma2 = sigma**2\n",
    "\n",
    "def metropolis(xold):\n",
    "    xtrial = np.random.random()    \n",
    "    xtrial = xold+(2*xtrial-1)*delta\n",
    "    weight = np.exp(-0.5*(xtrial**2-xold**2)/sigma2)\n",
    "    xnew = xold\n",
    "    if(weight >= 1): #Accept\n",
    "        xnew = xtrial\n",
    "    else:\n",
    "        r = np.random.random()\n",
    "        if(r <= weight): #Accept\n",
    "            xnew = xtrial\n",
    "    return xnew\n",
    "\n",
    "xwalker = 20.\n",
    "Nwarmup = 500000\n",
    "for i in range(Nwarmup):\n",
    "    xwalker = metropolis(xwalker)\n",
    "\n",
    "x[0] = xwalker\n",
    "tot = x[0]\n",
    "for i in range(1,N):\n",
    "    x0 = x[i-1]\n",
    "    for j in range(10):\n",
    "        x0 = metropolis(x0)\n",
    "    x[i] = metropolis(x0)\n",
    "\n",
    "binwidth=sigma/10\n",
    "pyplot.hist(x,bins=np.arange(-50,50., binwidth),density=True);\n",
    "\n",
    "norm = 1./(sigma*np.sqrt(2*np.pi))\n",
    "pyplot.plot(np.arange(-50.,50.,binwidth),norm*np.exp(-0.5*np.arange(-50.,50.,binwidth)**2/sigma2),ls='-',c='red',lw=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Challenge 9.2:\n",
    "\n",
    "- Modify the code above to study the equilibration \"time\" for different step size $\\delta$. \n",
    "\n",
    "- Analize the acceptance ratio in terms of $\\delta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}