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": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHvdJREFUeJzt3Xt0VfWd9/H3NwkJIDcJARHUhIsiSL2lqKPWC9ZBbaFW7GDVoX2cWh2trU7rg9NVnTqdjlqV6mhbXV5qO7aitI+lQr1iaxUvBLxgBCRclCBCuAhyDUm+zx97wzkcAzmQc7KTvT+vtc7K/v32PjnfzWZ9zs5v38zdERGRZCiIugAREWk7Cn0RkQRR6IuIJIhCX0QkQRT6IiIJotAXEUkQhb6ISIIo9EVEEkShLyKSIEVRF5CpT58+Xl5eHnUZIiIdypw5c9a4e1lLy7W70C8vL6eqqirqMkREOhQz+yCb5TS8IyKSIAp9EZEEUeiLiCSIQl9EJEEU+iIiCZJV6JvZGDNbaGY1ZjapmflfMLO5ZtZgZuMz5k00s0Xha2KuChcRkX3XYuibWSFwL3AOMBy4yMyGZyz2IfAN4HcZ7+0N3AScAIwCbjKzA1tftoiI7I9s9vRHATXuvsTd64HHgHHpC7j7Mnd/B2jKeO8/As+5+zp3Xw88B4zJQd27W7sWLr4Yjj8eRo7M+a8XEYmLbC7OGgAsT2vXEuy5Z6O59w7I8r3Z69YNpkyBxkYwgy1boGvXnH+MiEhH1y4O5JrZ5WZWZWZVdXV1+/4LSkqgoiKYdodFi3JboIhITGQT+iuAQ9LaA8O+bGT1Xne/390r3b2yrKzFW0c0b9iw1PTChfv3O0REYi6b0J8NDDWzCjMrBiYA07L8/c8AZ5vZgeEB3LPDvtxLD/0FC/LyESIiHV2Loe/uDcDVBGE9H3jc3avN7GYzGwtgZp83s1rgQuA+M6sO37sO+E+CL47ZwM1hX+4dcURqWqEvItKsrO6y6e4zgBkZfTemTc8mGLpp7r0PAQ+1osbsaE9fRKRF7eJAbk5kjuk3ZZ49KiIi8Qn9Pn2gtDSY3rIFVmR7rFlEJDniE/qgcX0RkRbEK/R12qaIyF7FN/S1py8i8hnxCn0N74iI7FW8Ql97+iIiexWv0K+ogE6dgukVK+DTT6OtR0SknYlX6HfqBEOGpNrvvx9dLSIi7VC8Qh94pqFnqqEhHhGR3cQu9BeXpu4Gcfe9f46wEhGR9id+od87dSfnQet0Va6ISLrYhf6S3qkHcw1eVxthJSIi7U/sQj99eGfQuhUMun4a5ZOmR1iRiEj7EbvQ39i5G3UH9AKgpHEHh2xYFXFFIiLtR+xCH+D9Pofumj58zYcRViIi0r7ENPQP2zU9VKEvIrJLLEN/YVroH1H3QYSViIi0L7EM/ffLUqF/+BqFvojITvEM/bQ9/cFraylqbIiwGhGR9iOWob+ppCu1PcoAKG5qoHz9RxFXJCLSPsQy9GH3vX2N64uIBOIb+hrXFxH5jPiGftq5+kco9EVEgBiHfvppm4dreEdEBIhx6NeUHkITBkD5+pWwdWvEFYmIRC+2ob+9UwnLDuwPQAEO8+dHXJGISPRiG/qw+8Fc3n03ukJERNqJWId++ri+Ql9EJOahn36uPtXV0RUiItJOxDr0F2p4R0RkN7EO/WUHHkx9QVHQ+PBD2LAh2oJERCKWVeib2RgzW2hmNWY2qZn5JWY2JZz/upmVh/2dzOwRM5tnZvPN7Ibclr93DYVFuz0+kXnz2vLjRUTanRZD38wKgXuBc4DhwEVmNjxjscuA9e4+BJgM3Br2XwiUuPtI4Hjg2zu/ENrKe30rUo233mrLjxYRaXey2dMfBdS4+xJ3rwceA8ZlLDMOeCScngqMNjMDHDjAzIqALkA9sDEnlWfpvb6DUo0332zLjxYRaXeyCf0BwPK0dm3Y1+wy7t4AbABKCb4ANgMrgQ+B2919XStr3ifv9Ru8a3re9L9RPml6W368iEi7ku8DuaOARuBgoAL4NzMblLmQmV1uZlVmVlVXV5fTAqr7pT7u8DUf6IEqIpJo2YT+CuCQtPbAsK/ZZcKhnJ7AWuDrwNPuvsPdVwOvAJWZH+Du97t7pbtXlpWV7fta7MXGzt2o7dEXgJLGBoau1YPSRSS5sgn92cBQM6sws2JgAjAtY5lpwMRwejww092dYEjnTAAzOwA4EViQi8L3Rfre/vBVS9v640VE2o0WQz8co78aeAaYDzzu7tVmdrOZjQ0XexAoNbMa4Dpg52md9wLdzKya4MvjYXd/J9cr0ZL0g7kjVi1u648XEWk3irJZyN1nADMy+m5Mm95GcHpm5vs2Ndff1qrTDuYOX70kwkpERKIV6ytyd3qvX+pc/eGrl4J7hNWIiEQnEaH/Ufcy1nfuDkCP7Zth2bJoCxIRiUgiQh+z3Q7m6iItEUmqZIQ+GVfm6nYMIpJQiQl97emLiCQo9LWnLyKSoNBfUjqQbUXFQaO2FtasibYgEZEIJCb0GwsKWVBWnuqYMyeyWkREopKY0Ad456Chqcbs2dEVIiISkUSF/tv9D0813ngjukJERCKSsNBP29N/4w1dmSsiiZOo0F/SewCfFncJGqtWBQd0RUQSJFGh31RQyDyN64tIgiUq9AHeyRziERFJkMSF/ltpB3Nn/f4vEVYiItL2Ehf66Xv6Iz9eBE1NEVYjItK2Ehf6H3Uvo+6AXgB0r98KCxdGXJGISNtJXOhjttsQjw7mikiSJC/0ybgyVwdzRSRBEhn6ujJXRJIqkaG/22mbb78N27dHV4yISBtKZOh/0qUHy3r1Dxr19UHwi4gkQCJDH+Ctg9OGeF59NbpCRETaUGJDv2rA8FTj5ZejK0REpA0lNvTnDDwy1XjlFd1xU0QSIbGhv7DPYWws7ho0Vq6EZcsirUdEpC0kNvSbCgp5c8CwVIeGeEQkARIb+gBVAzKGeEREYi7ZoT9QB3NFJFkSHfpv9T+CBgv/CaqrYf36aAsSEcmzRIf+1uLOVPcbnOrQ+foiEnOJDn2AOenj+hriEZGYS3zo7zaur4O5IhJzWYW+mY0xs4VmVmNmk5qZX2JmU8L5r5tZedq8z5nZq2ZWbWbzzKxz7spvvd3O4HnjjeBePCIiMdVi6JtZIXAvcA4wHLjIzIZnLHYZsN7dhwCTgVvD9xYB/wtc4e4jgNOBHTmrPgdWdy+FioqgsW0bzJ0bbUEiInmUzZ7+KKDG3Ze4ez3wGDAuY5lxwCPh9FRgtJkZcDbwjru/DeDua929MTel59App6Sm//a36OoQEcmzbEJ/ALA8rV0b9jW7jLs3ABuAUuBwwM3sGTOba2bXt77kPDjjjNT0iy9GV4eISJ7l+0BuEXAKcHH483wzG525kJldbmZVZlZVV1eX55I+65Q30m629ve/a1xfRGIrm9BfARyS1h4Y9jW7TDiO3xNYS/BXwUvuvsbdtwAzgOMyP8Dd73f3SnevLCsr2/e1aKXanv34sGe/oLFlC8yeTfmk6W1eh4hIvmUT+rOBoWZWYWbFwARgWsYy04CJ4fR4YKa7O/AMMNLMuoZfBqcB7+Wm9Nx69dDPpRozZ0ZXiIhIHrUY+uEY/dUEAT4feNzdq83sZjMbGy72IFBqZjXAdcCk8L3rgTsJvjjeAua6e7vchZ51WFroa1xfRGKqKJuF3H0GwdBMet+NadPbgAv38N7/JThts13bbU9/1ixKjvtOdMWIiORJ4q/I3Wl191IW9x4YNLZv57gVC6ItSEQkDxT6aV49dOSu6ZM+eDvCSkRE8kOhn2bWYUfvmj7pw3kRViIikh8K/TSvpe3pH7NyIWzeHGE1IiK5p9BPs65rTxgZBH+npsbgQi0RkRhR6Gc688zU9LPPRleHiEgeKPQzjRmTmn766ejqEBHJA4V+ptNOg87hLf/nz4cPPoi2HhGRHFLoZ+rSBU4/PdXW3r6IxIhCP0P5pOn8eEfa/eUU+iISIwr9Zvx1UGWq8fzzutWyiMSGQr8ZSw88OHWr5U2bYNasaAsSEckRhX5zzPjboONT7b/8JbpaRERySKG/B39ND32N64tITCj09+DVQz8HxcVB4513YEXmw8JERDoehf4ebCnuAqeemurQEI+IxIBCf2/OOSc1PS3zCZEiIh2PQn9vxo1LTT/7bHAmj4hIB6bQ34vyBxayoM9hQWP7dt2ATUQ6PIV+C54demKq8eST0RUiIpIDCv0WPHv4SanGU0/Bjh3RFSMi0koK/Ra8228wH3XvEzTWr9eDVUSkQ1Pot8RMQzwiEhsK/Sx8JvTdoytGRKQVFPpZeOOQo6BXr6CxfDm8+Wa0BYmI7CeFfhYaCovgS19KdfzhD9EVIyLSCgr9LH1786GpxpQpGuIRkQ5JoZ+lvw6qZGNx16CxeDFUVUVbkIjIflDoZ2l7UfHu5+w/9lh0xYiI7CeF/j7485FfSDWmTIGmpuiKERHZDwr9ffDKYUeztkuPoLFiBbz8crQFiYjsI4X+PmgoLOIvR5yc6tAQj4h0MFmFvpmNMbOFZlZjZpOamV9iZlPC+a+bWXnG/EPNbJOZfT83ZUdn2vDTdk2v/fWjDPnBnyKsRkRk37QY+mZWCNwLnAMMBy4ys+EZi10GrHf3IcBk4NaM+XcCsXj01OyBw1nZrRSA0q0b+YcP3o64IhGR7GWzpz8KqHH3Je5eDzwGjMtYZhzwSDg9FRhtZgZgZl8BlgLVuSk5Wm4FPHVk6jGK51e/GGE1IiL7JpvQHwAsT2vXhn3NLuPuDcAGoNTMugH/F/hx60ttP54cccau6XPenwUbNkRYjYhI9vJ9IPc/gMnuvtfnDJrZ5WZWZWZVdXV1eS6p9ar7Daa67yAAOjfU64CuiHQY2YT+CuCQtPbAsK/ZZcysCOgJrAVOAG4zs2XA94B/N7OrMz/A3e9390p3rywrK9vnlYjCEyPPSjUefji6QkRE9kE2oT8bGGpmFWZWDEwApmUsMw2YGE6PB2Z64FR3L3f3cuDnwE/d/Z4c1R6pJ0eczvbCoqDx+utQHYtDFiIScy2GfjhGfzXwDDAfeNzdq83sZjMbGy72IMEYfg1wHfCZ0zrj5pMuPXhuSNp99rW3LyIdgHk7u1tkZWWlV7XiZmblk6bnsJq9O23JHB554qag0bcv1NZCp05t9vkiIjuZ2Rx3r2xpOV2R2wp/Lz9m1zn7rF4NM2ZEW5CISAsU+q3QVFDIH0aOTnXcd190xYiIZEGh30pTPnc2BNehwdNPQ01NtAWJiOyFQr+Vlvc6CM49N2i4wy9/GW1BIiJ7odDPhavTLj146CHYvDm6WkRE9kKhnwtnnw1DhgTTn3wCjz4abT0iInug0M+FggK46qpU+5579OB0EWmXFPq58o1vQNfwwenz5vG1S26LtBwRkeYo9HOlVy+49NJdzf9TpYeriEj7o9DPgV1XAacd0D37/ddg4cKIKhIRaZ5CP5eOOoqZg4KroAtw+NnPIi5IRGR3Cv0c++WJ41ON3/4WPvooumJERDIo9HOkfNJ0yidNZ/bAEcw5eFjQWV8PP/95tIWJiKRR6OeaGb9K39v/1a+Cc/dFRNoBhX4ePD9kFItKw4eNffop/OIX0RYkIhJS6OeBWwH3nXDBrva6n9wKGzdGWJGISEChnydPDj8dKioA6L11I9x1V7QFiYig0M+bhsIiuOmmXe2NP7kF1q2LsCIREYV+fl18MYt7DwSgR/0WuOOOiAsSkaRT6OdTURGTT/l6qn3XXcFjFUVEIqLQz7Ppw05hfll50Ni8Gf77vyOtR0SSTaGfR+WTpuNWwJ2nXpLqvPdePVJRRCKj0G8Dzw05gdkDhgeNHTvgBz+ItiARSSyFflsw4z9H/0uq/eST8OKL0dUjIoml0G8j7/Q/fLf77XPttdDYGF1BIpJICv02dGK30WzpVBI03n47eIi6iEgbUui3oY979OFXJ6TdjO2GG2DNmugKEpHEUei3sftHnU9tj7KgsXYtfP/70RYkIomi0G9j2zp15sYvXpnqeOQRmDkzuoJEJFEU+hGYOWQUjE8b5rniCti2LbqCRCQxFPpRuftu6NkzmF60CH7yk2jrEZFEUOhHpX9/uOWWVPuWW2D27OjqEZFEUOhHpHzSdCqWDoBTTgk6Ghvhkktgy5ZoCxORWMsq9M1sjJktNLMaM5vUzPwSM5sSzn/dzMrD/i+a2Rwzmxf+PDO35XdsbgXwm9+wqbhL0PH++zxy6tconzQ92sJEJLZaDH0zKwTuBc4BhgMXmdnwjMUuA9a7+xBgMnBr2L8G+LK7jwQmAr/NVeGxUVHBj0dfvqs5ce50TlsyJ8KCRCTOstnTHwXUuPsSd68HHgPGZSwzDngknJ4KjDYzc/c33f2jsL8a6GJmJbkoPE6eGHkWzww9cVf79hmTYeXKCCsSkbjKJvQHAMvT2rVhX7PLuHsDsAEozVjmAmCuu2/P/AAzu9zMqsysqq6uLtva48OMG8Z8B/r1A6Bs8ycwYQI0NERcmIjETZscyDWzEQRDPt9ubr673+/ule5eWVZW1hYltTvruvbk66d/h0YLN8lLL8EPfwgEB301zi8iuZBN6K8ADklrDwz7ml3GzIqAnsDasD0Q+H/AP7v74tYWHGezyo/hzlMuTnXcdhv86U/RFSQisZNN6M8GhppZhZkVAxOAaRnLTCM4UAswHpjp7m5mvYDpwCR3fyVXRcdJ5h78L066kBcGfz7VcemlHF63rG2LEpHYajH0wzH6q4FngPnA4+5ebWY3m9nYcLEHgVIzqwGuA3ae1nk1MAS40czeCl99c74WMeJWwHXnXQfl5UHHp5/y0NSb6bN5vYZ4RKTVzN2jrmE3lZWVXlVVtd/vj0swDlu9lKmPXk+3+q0AvNn/CCZc9FMW3vHViCsTkfbIzOa4e2VLy+mK3HZqQd8Krh57PRQEm+jYlQu5Y/pkPW1LRFpFod+O/XXw52Hy5F3tLy18Ga66Ctxj8xeNiLQthX47V75iEA8f/+VUx333waRJ0M6G5USkY1Dot3dm3Dz6W/xxxBmpvttu419feyK6mkSkw1LodwBuBVx/znd5Nu1WDde/9Bvdg19E9plCv4NoKCziO2Ov5+XDjk51/uhHwVW7GuoRkSwp9DuQ7UXF/MsFP+Kl8mNTnT/9KVx7LTQ1RVeYiHQYCv0OZlunznzrgh/xfPpVu3fdBRddpOfsikiLFPod0PaiYq48/9+ZfsTJqc7HH4ezzoK1a6MrTETaPYV+B7WjsBPfGXv97qdzvvIKnHQSvPdedIWJSLum0O/AmgoKg6du3XknmAWdixbBqFEwdSqg2zKLyO4U+h2dGVx7LVeMm8SWTuFDyTZvhgsv5FcnXEBRox7EIiIpCv2YePqIk/nqJbezrFf/XX1XvPFHpj76Aw5b/9Fe3ikiSaLQj4GdwzcL+lYwduJkZg5K3WjvmJWLmPHwNfzg3O/qfH4RUejHzcbO3bhs/I389PRvUl9QBMABO7bxs7/cDeeeCx98EHGFIhIlhX4MuRVw/wkXcP6lt7O498DUjKefhhEj4H/+p9lbNOugr0j8KfRjrPqgIXxp4s95+Pgv00R4ds/mzXDNNXD88cHD10UkURT6Mbe1uDM/PuvbjL/kNhaVpj3f/u234bTTYMIEqKmJrkARaVMK/YSYO+BIzvvG3dx+6iVsLSpJzZgyBYYNg299i4M3ro6uQBFpE0VRFyBtp76oE/f8wwT+cNRoXl3/DDz2WDCjsREeeIAXC3/N748eA9ccS/ndbwKw7JbzIqxYRHJNe/oJtLJHGeWHXcJXL/4ZnHnmrv6Sxga+MfcpKC/njqfuYMSqxRFWKSL5oNBPsLkDj4QXXoAXXmDOwcNSM3bs4ILqF5n+6+/CGWfAn/+sB7KLxIR5O7tgp7Ky0quqqvb7/TrlcD+5c/qSKv71tScYVdvMDdsGDICJE+Gb34QhQ9q+PhHZKzOb4+6VLS2nPX0JmPHXwZ/naxffxrhL72DakV+gwdL+e6xYETywZejQ4KyfBx/UbZxFOiDt6cseHbxxNbO6L2DNPffRZ8uGzy5QWBgcExg/Hs4/H8rK2r5IEQGy39NX6EuLihobOHPxbC6c9xxfXDqn2fH9Rivgrf6Hc/wVX+fLNd1596DBLL31y838NhHJB4W+5MWy7x0Pv/sdc++4n+M+WrjH5dZ07Umfr5wHp54avI48Ego0miiSLwp9ybv+G+sY8/4sbtr2Hk2vzKKAvfxf6t0bTj45+AI4+WQ4+mg44IC2K1Yk5hT60maW3XIex17zO05d9ianLZnDaUvnNn8MIF1BQXAl8HHH7XqNfHIVn5YcoAvCRPaDQl8iY97E8FVLOHH5u3y+tprK2vda/hIIfdytNwedcGwwHDRsWOpn//6pR0KKyGco9KX9cGfwuloqa99jVG01FzSuhAUL9u2hLp07Q3l58KqoCF47pw89FPr00TEDSTSFvrRrXeq3cWTdUkasWsxRHy9m5KoaBq9dTsn+PtO3qAj69YODDw7+KujfPzVdVgalpcEXQ2lpcHyhSLedknjJNvSz+p9vZmOAu4BC4AF3vyVjfgnwG+B4YC3wT+6+LJx3A3AZ0Ahc4+7P7MN6SExtLe7M3AFHMnfAkbv6CpsaOeSTjxmytpYha5fveh2zrQ7Wr9/7L2xoCC4gW7EiuwJ69Qq+AHZ+GfTqBT16tPzq3h26doUuXaC4WENO0uG0GPpmVgjcC3wRqAVmm9k0d0+/Vv8yYL27DzGzCcCtwD+Z2XBgAjACOBh43swOd3fdyEU+o7GgkGW9B7Cs9wCeH3rCbvO6b9/MwA2rOOSTVQzcsJpDNny8q33QprX02rZp3z7sk0+C1+JW3FSuoCAI/51fAunTmX2dOwdfEjtfnTo1P723eUVFwQVxma899e/pVVCgL6sEy2ZPfxRQ4+5LAMzsMWAckB7644D/CKenAveYmYX9j7n7dmCpmdWEv+/V3JQvSfFpyQHM7zuI+X0HNTu/pKGesk3r6LtpPX03r6PfpnX0DdsHbt3AWX2LYO1a1n+4kp7bNu399NJsNTUFTyLbvLn1v6utFRR89svAbPfXzi+HfPbtrX+nzC+obOfty7Lt6TOuugrOy98ZbNmE/gBgeVq7FjhhT8u4e4OZbQBKw/7XMt47YL+rFdmD7UXF1PY6iNpeB7W4bEFTIz22b6b3lo302vopB27bSPftW+i2fQvd64Of3ep3trdmtLdQsqOezg31FDft5/GH9qCpKXjt2BF1JZJp7Ni8/vp2cTTLzC4HLg+bm8xsz5d6tqwPsKb1VXUoWudk0DonwZVX9uHKK/dnnQ/LZqFsQn8FkPZwVQaGfc0tU2tmRUBPggO62bwXd78fuD+bgltiZlXZHMGOE61zMmidkyHf65zNic2zgaFmVmFmxQQHZqdlLDMNmBhOjwdmenAu6DRggpmVmFkFMBR4Izeli4jIvmpxTz8co78aeIbglM2H3L3azG4Gqtx9GvAg8NvwQO06gi8GwuUeJzjo2wBcpTN3RESik9WYvrvPAGZk9N2YNr0NuHAP7/0v4L9aUeO+yskwUQejdU4GrXMy5HWd290VuSIikj+6WYmISILEJvTNbIyZLTSzGjObFHU9+WBmh5jZi2b2nplVm9l3w/7eZvacmS0Kfx4Yda25ZmaFZvammT0VtivM7PVwe08JTzKIFTPrZWZTzWyBmc03s5Pivq3N7Nrw//a7ZvZ7M+sct21tZg+Z2Wozezetr9ntaoG7w3V/x8yOa+3nxyL0024VcQ4wHLgovAVE3DQA/+buw4ETgavC9ZwEvODuQ4EXwnbcfBeYn9a+FZjs7kOA9QS3Aombu4Cn3X0YcDTB+sd2W5vZAOAaoNLdjyI4cWTnbV3itK1/DYzJ6NvTdj2H4KzHoQTXMv2ytR8ei9An7VYR7l4P7LxVRKy4+0p3nxtOf0oQAgMI1vWRcLFHgK9EU2F+mNlA4DzggbBtwJkEt/yAeK5zT+ALBGfG4e717v4JMd/WBCeXdAmv9+kKrCRm29rdXyI4yzHdnrbrOOA3HngN6GVm/Vvz+XEJ/eZuFRHr2z2YWTlwLPA60M/dV4azPgb6RVRWvvwcuB5oCtulwCfuvvM+CHHc3hVAHfBwOKz1gJkdQIy3tbuvAG4HPiQI+w3AHOK/rWHP2zXn2RaX0E8UM+sG/AH4nrtvTJ8XXhQXm1OyzOxLwGp3nxN1LW2sCDgO+KW7HwtsJmMoJ4bb+kCCPdsKgrvyHsBnh0FiL9/bNS6hn9XtHuLAzDoRBP6j7v7HsHvVzj/5wp+ro6ovD04GxprZMoJhuzMJxrp7hUMAEM/tXQvUuvvrYXsqwZdAnLf1WcBSd69z9x3AHwm2f9y3Nex5u+Y82+IS+tncKqLDC8eyHwTmu/udabPSb4MxEfhTW9eWL+5+g7sPdPdygu06090vBl4kuOUHxGydAdz9Y2C5mR0Rdo0muLI9ttuaYFjnRDPrGv5f37nOsd7WoT1t12nAP4dn8ZwIbEgbBto/7h6LF3Au8D6wGPhh1PXkaR1PIfiz7x3grfB1LsEY9wvAIuB5oHfUteZp/U8HngqnBxHcx6kGeAIoibq+PKzvMUBVuL2fBA6M+7YGfgwsAN4FfguUxG1bA78nOGaxg+Avusv2tF0BIzgzcTEwj+DMplZ9vq7IFRFJkLgM74iISBYU+iIiCaLQFxFJEIW+iEiCKPRFRBJEoS8ikiAKfRGRBFHoi4gkyP8H37Y1JZcneSYAAAAASUVORK5CYII=\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": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptance Ratio:  0.16499\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHuFJREFUeJzt3Xt0VeWd//H3NwnhLteA3CTc5FZEJaCOtVZoHawKnVZ/g9VVun6OtFbbTh1rcbrGVp126qVaO7W1/KqW2s6Il6lDlaJWtN4xARUJNwOiBFHD/Z4Q8v39sTecwyGQAznJTvb+vNY6K/t5zj453+3Gz9l59j7PNndHRESSIS/qAkREpPko9EVEEkShLyKSIAp9EZEEUeiLiCSIQl9EJEEU+iIiCaLQFxFJEIW+iEiCFERdQKaePXt6cXFx1GWIiLQqixYt2ujuRQ2t1+JCv7i4mLKysqjLEBFpVczs/WzW0/COiEiCKPRFRBJEoS8ikiAKfRGRBFHoi4gkSFahb2aTzWylmVWY2cx6nv+MmS02s1ozuyTjuelm9m74mJ6rwkVE5Ng1GPpmlg/cC1wAjAIuM7NRGat9AHwN+K+M13YHfgicAUwAfmhm3RpftoiIHI9sjvQnABXuvsbda4CHganpK7j7WndfAtRlvPbvgWfdfbO7bwGeBSbnoO5DbdoEl18O48bBmDE5//UiInGRzZez+gHr0tqVBEfu2ajvtf2yfG32OnWCOXNg/34wg927oUOHnL+NiEhr1yJO5JrZDDMrM7OyqqqqY/8FbdvCoEHBsju8+25uCxQRiYlsQn89MCCt3T/sy0ZWr3X3We5e4u4lRUUNTh1RvxEjUssrVx7f7xARiblsQr8UGGZmg8ysEJgGzM3y9z8NnG9m3cITuOeHfbmXHvorVjTJW4iItHYNhr671wLXEoT1cuARdy83s1vMbAqAmY03s0rgUuA3ZlYevnYzcCvBB0cpcEvYl3vDh6eWFfoiIvXKapZNd58HzMvouyltuZRg6Ka+1z4APNCIGrOjI30RkQa1iBO5OZE5pl+XefWoiIjEJ/R79oQePYLl3bthfbbnmkVEkiM+oQ8a1xcRaUC8Ql+XbYqIHFV8Q19H+iIih4lX6Gt4R0TkqOIV+jrSFxE5qniF/qBB1OSFXz1Yvx527Ii2HhGRFiZeod+mDe9365Nqr1oVXS0iIi1QvEIfWNM9beZmDfGIiBwidqG/ukfabBAKfRGRQ8Qv9LunzeSsa/VFRA4Ru9DX8I6IyJHFLvQPGd5ZtSq4haKIiAAxDP3t7TpB795Bo7oa3nsv2oJERFqQ2IU+AKNHp5bLy6OrQ0SkhYll6D+4rWOqodAXETkolqG/sufAVGPp0ugKERFpYWIZ+quKFPoiIvWJZ+inH+mvWAH79kVXjIhICxLL0N/ZtgOVJxQFjX374N13oy1IRKSFiGXoQ8bRvoZ4RESAOIe+xvVFRA4T39DveVKqodAXEQFiHPq6bFNE5HCxDf2KHgOow8JGBezZE21BIiItQGxDv7pNW9YeuIuWOxdePSvagkREWoDYhj4cejJ3+Mb3I6xERKRliHXop4/rn1yl0BcRiXXop1+rf/LGDyKsRESkZYh16K8sSg99HemLiMQ69Nd260tNXgEA/bdXwbZtEVckIhKtrELfzCab2UozqzCzmfU839bM5oTPLzSz4rC/jZnNNrN3zGy5md2Y2/KPrja/4NDbJ77zTnO+vYhIi9Ng6JtZPnAvcAEwCrjMzEZlrHYlsMXdhwJ3A7eF/ZcCbd19DDAO+PqBD4TmsqzXoFTjrbea861FRFqcbI70JwAV7r7G3WuAh4GpGetMBWaHy48Bk8zMAAc6mlkB0B6oAbbnpPIsLes1ONV4883mfGsRkRYnm9DvB6xLa1eGffWu4+61wDagB8EHwC5gA/ABcKe7b25kzcdkWe8hqYaO9EUk4Zr6RO4EYD/QFxgE/IuZDc5cycxmmFmZmZVVVVXltIDy3mlvt3SpbqgiIomWTeivBwaktfuHffWuEw7ldAE2AV8B5rv7Pnf/BHgFKMl8A3ef5e4l7l5SVFR07FtxFNvbdaLyhF5Bo6YGli3L6e8XEWlNsgn9UmCYmQ0ys0JgGjA3Y525wPRw+RJggbs7wZDORAAz6wicCazIReHH4pCjfQ3xiEiCNRj64Rj9tcDTwHLgEXcvN7NbzGxKuNr9QA8zqwCuAw5c1nkv0MnMygk+PB509yW53oiG6GSuiEigIJuV3H0eMC+j76a05b0El2dmvm5nff3NrVwnc0VEgJh/I/eAZb0zrtV3j64YEZEIJSL0P+xcxJZ2nYPGtm2wdm2k9YiIRCURoY/ZoSdzNa4vIgmVjNAn42SuxvVFJKESE/o60hcRSVDo60hfRCRBob+mR39o1y5oVFbCxo3RFiQiEoHEhP7+vHw45ZRUx6JF0RUjIhKRxIQ+ACVp0/6UlkZXh4hIRJIV+hMmpJbfeCO6OkREIpKs0B8/PrX8xhv6Zq6IJE6iQn/w7yrYUdg+aHz8cXBCV0QkQRIV+nV5+bxz4rBUh8b1RSRhEhX6AEv6pIW+xvVFJGESF/pv9Tk51VDoi0jCJC70DznSLyuDurroihERaWaJC/0POxdR1bFr0Nixg0kzZkVbkIhIM0pc6GN2yBDP2I9WRViMiEjzSl7oA0vSruAZu0GhLyLJkcjQfzv9SF+hLyIJksjQTz+ZO/KT96C6OsJqRESaTyJDf2v7E1jbtQ8AbffXwttvR1yRiEjzSGToA7zVN+16/ddei64QEZFmlNjQL+s3KtV4+eXoChERaUaJDf1F/UemGq+8ohk3RSQREhv6K3sOZHthh6CxYQOsXRtpPSIizSGxoV+Xl8+b/UakOjTEIyIJkNjQByjrlzHEIyISc8kO/f46mSsiyZLo0H+rz3BqLfxPUF4OW7ZEW5CISBNLdOjvKWxHee8hqQ5dry8iMZfo0AdYlD6uryEeEYm5xIf+IeP6OpkrIjGXVeib2WQzW2lmFWY2s57n25rZnPD5hWZWnPbcKWb2mpmVm9k7ZtYud+U3XvoVPHtffR1qaiKsRkSkaTUY+maWD9wLXACMAi4zs1EZq10JbHH3ocDdwG3hawuAPwDfcPfRwGeBfTmrPgc+6dwDBg0CoF1tDSxeHHFFIiJNJ5sj/QlAhbuvcfca4GFgasY6U4HZ4fJjwCQzM+B8YIm7vw3g7pvcfX9uSs+hT386tfy3v0VXh4hIE8sm9PsB69LalWFfveu4ey2wDegBnAy4mT1tZovN7IbGl9wEzjsvtfz889HVISLSxAqa4fd/GhgP7AaeM7NF7v5c+kpmNgOYAXDSSSc1cUn1SA/9l14KxvULC5u/DhGRJpbNkf56YEBau3/YV+864Th+F2ATwV8FL7r7RnffDcwDTs98A3ef5e4l7l5SVFR07FvRWMXFB8f12b0bSkubvwYRkWaQTeiXAsPMbJCZFQLTgLkZ68wFpofLlwAL3N2Bp4ExZtYh/DA4F1iWm9Jzp3jmU8zpNDTVsWBBdMWIiDShBkM/HKO/liDAlwOPuHu5md1iZlPC1e4HephZBXAdMDN87RbgLoIPjreAxe7+VO43o/FeHXhKqqFxfRGJqazG9N19HsHQTHrfTWnLe4FLj/DaPxBcttmivXZSWui/+irs3QvtWtRXCkREGi3x38g94JPOPVjdvX/QqK7WPDwiEksK/TSvnTQm1dC4vojEkEI/zasDx6YaGtcXkRhS6Kd5Pf1If+FC2LUrumJERJqAQj/N5g5dYEwY/LW1wRe1RERiRKGfaeLE1PIzz0RXh4hIE1DoZ5j+YfdUY/786AoREWkCCv0Mrw/4FHsLwnl3li+H99+PtiARkRxS6GeobtOW1wekndDV0b6IxIhCvx5/G5w2J5xCX0RiRKFfjxcGlxxc3jnvad1CUURiQ6Ffj/e69eWDLr0B6FSzJ5iLR0QkBhT69THjb4PHpdp/+Ut0tYiI5JBC/wheSA99jeuLSEwo9I/gtZNOoTo/nHl6yRJYn3mzMBGR1kehfwS7C9tT2n90qkNDPCISAwr9ozhkiGdu5h0iRURaH4X+UTw77MxU45lnYOfO6IoREckBhf5RvN+tLyt6Dgwa1dWagE1EWj2FfgOeST/af+KJ6AoREckBhX4Dnjn5rFTjySdh377oihERaSSFfgOW9h7Ch517Bo0tW3RjFRFp1RT6DTHTEI+IxIZCPwuHhb57dMWIiDSCQj8Lbwz4FHTtGjTWrYM334y2IBGR46TQz0JtfgFcdFGq4/HHoytGRKQRFPpZ+vqukw4ur/3VgxriEZFWSaGfpRcGl7C9sAMAxVs3QFlZxBWJiBw7hX6WqgsKD71m/+GHoytGROQ4KfSPwZ9HfibVmDMH6uqiK0ZE5Dgo9I/BKwPHsqn9CUFj/Xp4+eVoCxIROUYK/WNQm1/AX4afnerQEI+ItDJZhb6ZTTazlWZWYWYz63m+rZnNCZ9faGbFGc+fZGY7zez63JQdnbmjzk01Hn0U9u2jeOZT0RUkInIMGgx9M8sH7gUuAEYBl5nZqIzVrgS2uPtQ4G7gtozn7wJiceup0v6j2NCpR9DYuBEWLIi2IBGRY5DNkf4EoMLd17h7DfAwMDVjnanA7HD5MWCSmRmAmX0ReA8oz03J0XLL48mR5xxs/+m6n0ZYjYjIsckm9PsB69LalWFfveu4ey2wDehhZp2A7wM3N77UluOJ0ecdXL5g1at0rt4VYTUiItlr6hO5PwLudvej3mfQzGaYWZmZlVVVVTVxSY1X3nsI5b0GA9CutoaLl78YcUUiItnJJvTXAwPS2v3DvnrXMbMCoAuwCTgDuN3M1gL/DPyrmV2b+QbuPsvdS9y9pKio6Jg3IgqPjvncweVLl/w1wkpERLKXTeiXAsPMbJCZFQLTgLkZ68wFpofLlwALPHCOuxe7ezHwc+An7v7LHNUeqSdGf5bq/AIATtuwEspjccpCRGKuwdAPx+ivBZ4GlgOPuHu5md1iZlPC1e4nGMOvAK4DDrusM262tj+BZ4emzbP/4IPRFSMikiXzFjZbZElJiZc1YjKz5rxm/tw1i5j96A+DRq9eUFkJbdo02/uLiBxgZovcvaSh9fSN3EZ4qfjU1DX7n3wC8+bpi1oi0qIp9BuhLi+fx8dMSnX85jfRFSMikgWFfiPNOeV86rCgMX8+A7d8GG1BIiJHodBvpHVdT+T5IeEwmjtXvDkv2oJERI5CoZ8Dvz89df/c/7PkWdilb+iKSMuk0M+BFwedxnvd+gDQpXoX/PGPEVckIlI/hX4OuOXx0Gmpo31++UvdOF1EWiSFfo48NmYSu9u0DRrvvAMvvRRtQSIi9VDo58j2dp34U9rsm/z859EVIyJyBAr9HJqddkKXJ56AlSujK0ZEpB4K/RxaVVTMgsGpyze5445oCxIRyaDQz7Ffn3lJqvHQQ/ChvqwlIi2HQj/HSvuPhrPOCho1NRrbF5EWRaGfa2bw/e+n2vfdB1u3RlePiEgahX5TuPhiGDkyWN6xA371q2jrEREJKfSbQl4e3HBDqn3XXbB9e3T1iIiEFPpNoHjmUwwt78YHXXoHHZs2wT33RFuUiAgK/SZTm1/APWd/JdVx552weXN0BYmIoNBvUk+M/iwMHx40tm+Hn/0s0npERBT6TWh/Xj7cfHOq4557gtsqiohERKHf1C69FMaMCZZ37YL/+I9o6xGRRFPoN7W8PLj11lT73nuhoiK6ekQk0RT6zWHKFEr7jQqW9+2D730v2npEJLEU+s3BjFsn/VOq/cQT8Pzz0dUjIoml0G8mS/qczOPp8+1/97uwf390BYlIIin0m9Edn5meurvW22/DAw9EW5CIJI5Cv4kVz3yK4plPAfDRCT2574y0qZdvvBE2boyoMhFJIoV+M5s14R+oPKEoaGzaBNdfH21BIpIoCv1mtrdNO276/NWpjtmzYcGC6AoSkURR6EdgwdAJPDX87FTHN74Be/dGV5CIJIZCPyI/+tzX2d62Y9B4913493+PtiARSQSFfkSqOnXntnOnpzp++lMoLY2uIBFJBIV+hP7r1Mm80T/8pu7+/XDFFbB7d7RFiUisZRX6ZjbZzFaaWYWZzazn+bZmNid8fqGZFYf9nzezRWb2TvhzYm7Lb93c8rjuwuvYWdg+6Fi16tA7bomI5FiDoW9m+cC9wAXAKOAyMxuVsdqVwBZ3HwrcDdwW9m8ELnb3McB04KFcFR4XlV1P5OZJM1Id994L8+dHV5CIxFo2R/oTgAp3X+PuNcDDwNSMdaYCs8Plx4BJZmbu/qa7fxj2lwPtzaxtLgqPk0fHfI6nh52Z6vja12DDhsjqEZH4yib0+wHr0tqVYV+967h7LbAN6JGxzpeBxe5enfkGZjbDzMrMrKyqqirb2uPDjBsnf4uqjl2D9scfw7RpUFsbbV0iEjvNciLXzEYTDPl8vb7n3X2Wu5e4e0lRUVFzlNTibO7Qhe9cdH0w/z7Aiy/CD34QbVEiEjvZhP56YEBau3/YV+86ZlYAdAE2he3+wJ+Ar7r76sYWHGevFp/KHWdfnuq4/Xau+vK/RVeQiMRONqFfCgwzs0FmVghMA+ZmrDOX4EQtwCXAAnd3M+sKPAXMdPdXclV0nP3qrEt5bsj4g+27nvwZLF0aYUUiEicNhn44Rn8t8DSwHHjE3cvN7BYzmxKudj/Qw8wqgOuAA5d1XgsMBW4ys7fCR6+cb0WMHLiMc12X3gB0rtkDF10UjPOLiDSSuXvUNRyipKTEy8rKjvv1B6Yxbu1GfPIej/3xBjrV7Ak6zjgjuNtW+/bRFiYiLZKZLXL3kobW0zdyW6gVvQZx7ZQb2G/hLlq4kCdPO1932xKRRlHot2AvDBnPrRNT99a9aOXLcM010ML+OhOR1kOh38L9btzFPDju4lTHb34DM2cq+EXkuBREXYA0wIxbJl1Fl707+VL580Hf7bdDly4Ubx8LwNqfXhhhgSLSmuhIvxVwy+OGC77DM+lTNfzgB1z76sPRFSUirZJCv5WozS/gW1NugEmTDvZd/9IfuP7F32uoR0SyptBvRaoLChkx9mpeLD7tYN+1rz3CA+O/CHV1EVYmIq2FQr+V2dumHVd9+d/4a9q3dv/vorlw2WW6z66INEih3wpVFxRy9T/866E3V3/kEUqHjoNNm6IrTERaPIV+K7Uvvw3fmnLDIZdzjl+/jDXDToFlyyKsTERaMoV+K1aXl8/Nk2Zw68R/og4DYPCWD2HCBHjssYirE5GWSKHf2plx//gv8s0vzmR3m/CmZLt2waWXBvfb3bcv2vpEpEVR6MfE/OFn86Ur7mRt1z6pzjvugLPPhoqK6AoTkRZFoR8jK3oNYsr0u+ELX0h1lpaya9QYeOABXc8vIgr9uNnerhP8+c/BVA1t2gDQcd9euPLK4MPg/fcjrlBEoqTQj6O8PPje9+D111ndvX+qf/58dg0bAf/5n5qiWSShFPoxVvzIBi6a/vPgsk4Lru7puG8vfPvbMG5ccPN1EUkUhX7M7Slsx82f+zpfuvx23u2Rdn/7t9+Gc8+FadN0olckQRT6MVQ886nDbhu5uN9ILvzaL7jznCvYU9A29cScOTBiBFx1FXzwwcHXi0g8KfQTpKagDb/8u2lMvOq+4Aj/gP374be/pXrwEH437mJ67dBUDiJxpdBPoA0nFFE88Aq+dPkdvDLwlIP9bffX8rXFT/LyfVfCV78Kb74ZYZUi0hQU+gm2uP9ILp/2Ey6b9mMW9R1xsL+wrhYeeghOPx3OOy+4BFRX+4jEgm6XKLw2cCxfvuIUPrumjG++/igTKtMmbHvhheDRrx9Mn865G4t5v1tfQLdpFGmNFPoSMOOFIeN5Ych4xn64kv/dVwqPPpo6wl+/Hn7yE/4GLBzwKR4fPRE2nQk9ekRatogcGw3vyGHe7juc4oFX8Hcz/h/3nfFlNnbocsjzZ6xbyu3zfwG9e8P558OsWVBVFVG1InIszFvYfCwlJSVeVlZ23K/X5Ya5V7C/lomrS5lVu4TaJ5+iwA+/NeN+y+OtPicz7htfgQsuCM4H5OmYQqS5mNkidy9paD39XykNqs0v4JmTz4K5cznrm7O59bwrWdx3+CHr5Hsd4z5cATfdBOPHw4knwuWXw333QXm57uEr0kLoSF+OW5/tVUxe9So/3LuMuldeJY+j/Fvq3j2Y5vmcc4KfY8dCx47NV6xIzGV7pK/Ql5zotnsb56x9k190/ZiNj8+l5+5tR39BXl7wTeDTT089Tj0VunQ5+utEpF4KfYmMeR2jPl7DmeuWMr6ynJLKZQ1/CBzQty+MHBk8RoxI/ezT5+CkcSJyOIW+tBzuDNlcSUnlMiZUlnPKhncZsqny6MNBmdq1g+Li4DFoUPA4sHzSSdCzp04cS6JlG/q6Tl+anhmrewxgdY8BzBn79wC0r9nLyKr3GP3xaj710WrGfFzBqK3roaam/t+xdy+sWBE86lNQEFxC2rdv8FdBnz6p5aKi4PsEPXsGP7t3D9YXSaCs/uWb2WTgHiAf+K27/zTj+bbA74FxwCbgH919bfjcjcCVwH7g2+7+dM6ql1ZrT2E7FvcbyeJ+Iw/25dftZ8DWjxi6qZKhm9YdfBRv+ZCue3ce/RfW1gZfIFu/PrsCunYNPgAOfBh07QonnNDwo3Nn6NAB2reHwkINOUmr02Dom1k+cC/weaASKDWzue6e9l19rgS2uPtQM5sG3Ab8o5mNAqYBo4G+wF/N7GR310Qucpj9efms7d6Ptd378ddhZxzyXOfqXfTf9jEDtn5M/22fMGDbRwfbI/Zvhy1bju3Ntm4NHqtXH3+9lkd+h/apD4H2acuZfe3aBR8SBx5t2tS/fLTnCgogP//wx5H6j/TIy9OHVYJlc6Q/Aahw9zUAZvYwMBVID/2pwI/C5ceAX5qZhf0Pu3s18J6ZVYS/77XclC9JsaNtR5b3GszyXoPrfb5tbQ1FOzfTa+cWeu3aTO+dm+kVtrvt2Ua3PTvouncHQ9hD3eYtx3Y+4QjyvQ527QoerU1e3uEfBmaHPg58ODRl39H6D8j8gMr2uWNZtyW9xzXXwIVNN69VNqHfD1iX1q4EzjjSOu5ea2bbgB5h/+sZr+133NWKHEF1QSGVXU+ksuuJDa6bV7efE6p30X33drru2UG3vdvpXL2bTtW76VwT/OxUc6C9J6O9m7b7amhXWxPMRtpa1dUFj337oq5EMk2Z0qS/vkWczTKzGcCMsLnTzFY24tf1BDY2vqpWRducDNrmJLj66p5cffXxbPPAbFbKJvTXA2k3V6V/2FffOpVmVgB0ITihm81rcfdZwKxsCm6ImZVlc9lSnGibk0HbnAxNvc3ZXNhcCgwzs0FmVkhwYnZuxjpzgenh8iXAAg++ADAXmGZmbc1sEDAMeCM3pYuIyLFq8Eg/HKO/Fnia4JLNB9y93MxuAcrcfS5wP/BQeKJ2M8EHA+F6jxCc9K0FrtGVOyIi0clqTN/d5wHzMvpuSlveC1x6hNf+GPhxI2o8VjkZJmpltM3JoG1Ohibd5hY3DYOIiDQdTVYiIpIgsQl9M5tsZivNrMLMZkZdT1MwswFm9ryZLTOzcjP7Ttjf3cyeNbN3w5/doq4118ws38zeNLMnw/YgM1sY7u854UUGsWJmXc3sMTNbYWbLzeysuO9rM/tu+G97qZn9t5m1i9u+NrMHzOwTM1ua1lfvfrXAL8JtX2Jmpzf2/WMR+mlTRVwAjAIuC6eAiJta4F/cfRRwJnBNuJ0zgefcfRjwXNiOm+8Ay9PatwF3u/tQYAvBVCBxcw8w391HAGMJtj+2+9rM+gHfBkrc/VMEF44cmNYlTvv6d8DkjL4j7dcLCK56HEbwXaZfN/bNYxH6pE0V4e41wIGpImLF3Te4++JweQdBCPQj2NbZ4WqzgS9GU2HTMLP+wIXAb8O2ARMJpvyAeG5zF+AzBFfG4e417r6VmO9rgotL2off9+kAbCBm+9rdXyS4yjHdkfbrVOD3Hngd6GpmfRrz/nEJ/fqmioj1dA9mVgycBiwEerv7hvCpj4DeEZXVVH4O3AAcuNFuD2Crux+YByGO+3sQUAU8GA5r/dbMOhLjfe3u64E7gQ8Iwn4bsIj472s48n7NebbFJfQTxcw6AY8D/+zu29OfC78UF5tLsszsIuATd18UdS3NrAA4Hfi1u58G7CJjKCeG+7obwZHtIIJZeTty+DBI7DX1fo1L6Gc13UMcmFkbgsD/o7v/T9j98YE/+cKfn0RVXxM4G5hiZmsJhu0mEox1dw2HACCe+7sSqHT3hWH7MYIPgTjv688B77l7lbvvA/6HYP/HfV/DkfdrzrMtLqGfzVQRrV44ln0/sNzd70p7Kn0ajOnA/zZ3bU3F3W909/7uXkywXxe4++XA8wRTfkDMthnA3T8C1pnZ8LBrEsE322O7rwmGdc40sw7hv/UD2xzrfR060n6dC3w1vIrnTGBb2jDQ8XH3WDyALwCrgNXAD6Kup4m28dMEf/YtAd4KH18gGON+DngX+CvQPepam2j7Pws8GS4PJpjHqQJ4FGgbdX1NsL2nAmXh/n4C6Bb3fQ3cDKwAlgIPAW3jtq+B/yY4Z7GP4C+6K4+0XwEjuDJxNfAOwZVNjXp/fSNXRCRB4jK8IyIiWVDoi4gkiEJfRCRBFPoiIgmi0BcRSRCFvohIgij0RUQSRKEvIpIg/x8OJgyqWdS2VgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 100000\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": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD8CAYAAABpcuN4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XucTdX7wPHPMzNm3MmghBqhCyVqQqRCQnL5RoUulPKrSIlqVFRUqNyTUgrfb4V8uyiVfNG9ZEiEZBJyKfdxnzEzz++PvefMmWkuBzOzz5x53q/XvGavtdfenl2c56y9115LVBVjjDEmJ2FeB2CMMSa4WaIwxhiTK0sUxhhjcmWJwhhjTK4sURhjjMmVJQpjjDG5skRhjDEmV5YojDHG5MoShTHGmFxFeB1AfqhcubLGxMR4HYYxxhQpy5cv362qVfJqFxKJIiYmhvj4eK/DMMaYIkVENgfSzm49GWOMyZUlCmOMMbmyRGGMMSZXliiMMcbkyhKFMcaYXFmiMMYYkytLFMYYY3JlicIYY0yuAkoUItJORNaLSIKIxGWzP0pEZrv7l4pIjFvfRkSWi8hq93crv2MudesTRGSiiIhbX0lEForIBvf3aflzqcYYY05Gnm9mi0g4MBloA2wFlonIPFVd69esD7BPVeuISHdgNHAzsBvoqKrbReRCYAFQ3T1mCnA3sBT4BGgHfArEAYtUdZSblOKAR0/9Uo3JXzFx87Ot3zSqQyFHYkzBCmQKj8ZAgqpuBBCRWUBnwD9RdAaecrfnAi+JiKjqT35t1gClRCQKqASUV9Uf3HPOBLrgJIrOwNXuMTOAL7BEYYqq5GTYsAHWrYM//oAzz4R69eC886B0aa+jMyYggSSK6sCffuWtQJOc2qhqiogkAtE4PYp0XYEVqpokItXd8/ifM72ncbqq7nC3/wJOzy4oEekL9AU466yzArgMYwpegx2/wdAfYO1a52fDBkhN/WdDEYiJcZJGvXrQujVce61Tb0yQKZRJAUWkPs7tqGtP5DhVVRHRHPZNBaYCxMbGZtvGmMLSZMtq+n83mxabVwZ2gKrTw/jjD5g/H154ARo1gieegC5dIMzGmZjgEUii2AbU9CvXcOuya7NVRCKACsAeABGpAbwP3K6qv/u1r5HDOf8WkWqqukNEqgE7T+B6jCk8qrTY9BP9v5tNk61rcm539tlOr+Gcc2DbNqenkZAAaWmZ2/30E3Tt6rR9/HG46SaICIkJnk0RF8jfwmVAXRGphfNh3h3omaXNPKAX8D3QDVjs9gYqAvOBOFX9Nr2xmwQOiEhTnIfZtwOTspxrlPv7w5O9OGMKzGef8cG/B9Fwx2+ZqlMkjIiePaBt24xnEWXL/vP4pCT47TcnaXz9NbzxBhw96uxbuxZuuQWefNLpYdx+u92SMp4S1bzv2ojIdcB4IBx4Q1WfFZHhQLyqzhORksC/gUbAXqC7qm4UkSeAIcAGv9Ndq6o7RSQWmA6UwnmIfb+bXKKBOcBZwGbgJlXdm1t8sbGxautRmEKRkgKPPgpjx2aqTg6L4L8XtmJK0xvZclq1HA/PcUTUzp3OOSdPhkOHMu+74QaYPh3KlTvF4I3JTESWq2psnu0CSRTBzhKFKRS7dsHNN8OSJb6qpPASzLr4Wl5t0pXt5avmeYo8h87u3QsTJ8KECbB/f0b9BRfA++87PRRj8kmgicKemBmTh5i4+XTsNZ5ttetnShIL6zThinum8WSbewNKEgGpVAmeego2bYL+/TPq162Dxo1h3rz8+XOMOQGWKIzJw42rFjL3rUeofnAXAGkIL7a4lb43PM6uspUK5g+tUAEmTYKZM6FkSafuwAHo3Nl5dpH1QbgxBchuPRmTk+PH4cEH4eWXfVUHosowoONgvqh9Wb7+UbneklqxwnlOsdlveeMOHeCtt5yEYsxJsltPxpyKtDTo1StTkvi18tl07DUu35NEni65BOLjnZfy0s2fD+3bw+HDhRuLKZYsURiTlSrcfz+8846v6uPzW3DDbS+y+bQzvYmpcmX47DN4+OGMuu+/d3oaycnexGSKDbv1ZExWw4bBiBG+4r8bXcfQNvd68i5DtrekJk6EBx7IKN90E7z9NoSHF15gJiTYrSdjTsb48ZmSBD16MKzNPcH1wtuAAc4D7XRz5kC/fk5PyJgCYInCmHQzZsDAgRnl9u1hxgxUgvCfyZNPZh4+++qrzrQfxhQAm0jGGIAPP4Q+fTLKzZvD3LlQooR3MZHHmhcTJsC+fc7oJ4CRIyE6GgYNKsQITXEQhF+VjClkX33lvHGdPh34xRfDxx8H/3oRYWHw5ptw/fUZdYMHO9N9GJOPLFGY4m37dujWzZmkD6B2bWd0UcWK3sYVqBIlnGcULVpk1P3f/znvXhiTTyxRmOIrJQV69nTmcAI4/XRYuBDOOMPbuE5UqVLw0UfQoIFTTk52RkIlJnoblwkZlihM8TV8OHz5pbMdFua8N1GrlrcxnawKFZxnKukzzP7+O9x9t42EMvnCEoUpnhYuhGeeySg/+SS0bOldPPmhbl147bWM8rvvwpQp3sVjQoaNejLFz/btzsJA6d+2W7eGxx/PcYRRkXLzzU4vKT1BDBwITZs604AYc5IsUZjiJetziTPOcIaXFrG3mnNLalGl2/Be1QXU37kRkpPZ3LoD1/eewOpxNxVihCaUBHTrSUTaich6EUkQkbhs9keJyGx3/1IRiXHro0VkiYgcEpGX/NqXE5GVfj+7RWS8u6+3iOzy23dX/lyqMWR6LpEqYfS4egAx4+JDozfhSoqIpF/nRzkYWQqAs/f/xahPJ9rzCnPS8kwUIhIOTAbaA/WAHiJSL0uzPsA+Va0DjANGu/XHgKHAYP/GqnpQVRum/+AsefqeX5PZfvtfP5kLM+YfsjyXmNC8B9+f3cDDgArOpkrViWs3wFfusP7bTDPhGnMiAulRNAYSVHWjqiYDs4DOWdp0Bma423OB1iIiqnpYVb/BSRjZEpFzgarA1yccvTGB2rcPbr/d963667Mb8tLloX0rZv4FLZjZyG9SwYceclbKM+YEBZIoqgN/+pW3unXZtlHVFCARiA4whu44PQj/fnFXEVklInNFpGaA5zEmZ4MHw19/AbCrTEUGdhxEWljRei5xMp5t1YdfTq/tFJKTnSGztjqeOUHBMDy2O/COX/kjIEZVGwALyeipZCIifUUkXkTid6U/mDQmO4sWwRtv+IqPX9uP3WVO8zCgwpMUEcng6x7keHpS/PZbGzJrTlggiWIb4P+tvoZbl20bEYkAKgB78jqxiFwMRKjq8vQ6Vd2jqu58CrwOXJrdsao6VVVjVTW2SpUqAVyGKZYOH3a+Rafr1o3Pz73cu3g88GvVWkxp0i2jIi4OtmzxLiBT5ASSKJYBdUWklohE4vQA5mVpMw/o5W53AxZrYCsi9SBzbwIRqeZX7ATYTVVz8oYNgz/+cLZPOw0mTfI2Ho+81Kw7nH++Uzh0CO65x0ZBmYDlmSjcZw79gQU4H9pzVHWNiAwXkU5us2lAtIgkAA8BviG0IrIJGAv0FpGtWUZM3USWRAEMEJE1IvIzMADofVJXZsyyZc5CROnGji168zjlk+SIEjBtWsYCTJ9+6qyKZ0wAbClUE5LqPvwB82YM5IJdmwD4KqYRt980PLhWqitkm0Z1cNYCf8l9pSk62hkFZbduiy1bCtUUa/f8MNeXJI6UiOKxtv2KdZIA523u+hFXsq2cmxj27OGDq2/0NihTJFiiMKFn3Tr6fz/bVxzT4ja2Viyet5yyOhxVmsfb9vOVu6z9EuaHzlvppmBYojChJS0N7rqLqNQUAFZWO5c3L+3ocVDB5YvasbxX32+m3HvvhYMHvQvIBD1LFCa0zJgB330HQHJYBI+0H1AsXqw7USNa3cWeUuWdwp9/Zp5y3ZgsLFGY0HHgAAwZ4itObXIDv1WJ8S6eILavdAVGtPZ7v2TcONiwwbuATFCzRGFCxzPPwN9/A7CjbDSTm4b2XE6n6oN6V7P8TPfdiuPHYdAgbwMyQcvWozBFQl7TgMfs3cbn08YR6ZZHtryDo5ElCz6wokyEp675Pz6aOdApf/QRLFgAbdt6G5cJOtajMCHhicWvE5nmPMCOr34B8y64yuOIiobV1erCHXdkVAwc6PQujPFjicIUeVdtXM41vy8DIA3nW3Jxf2fihDz3HJQr52yvW2frVph/sERhirSI1BSGLnrNV373omv45Yw6HkZUBJ1xBgwdmlF+8smMpWKNwRKFKeJuXzGfOnu3AnAwshQvXnm7xxEVUQMGQB03wSYmZk4cptizRGGKrEpHEnnw24yJ7SY268GussVjnYn8FBM3n5gn/8edF/f01aW9OhVWrvQwKhNMLFGYImvwV/+mfNJhADaedibTY+0N7FOxuPZlfFnrEgDCUHjgAZuK3ACWKEwRdcHOjXT/eYGvPKL13RwPL+FhRCFAhOGt7s5YDe+rr2DuXG9jMkHBEoUpkuK+mO586wW+qHUpS87Jc6ZkE4DfK9dk5iXXZ1QMGWLDZY0lClP0NNu0kqv+WAFAqoTxbMs7bThsPprQvAeJUWWcwu+/w2uv5X6ACXmWKEyRIppG3JfTfeW5F7ZmQ5WzvQsoBB0oWZbJl/tNf/L0087yqabYCihRiEg7EVkvIgkiEpfN/igRme3uXyoiMW59tIgsEZFDIvJSlmO+cM+50v2pmtu5jAG47tdvafBXAgDHIiIZf0XPPI4wJ2PGJddDjRpOYedOZxlZU2zlOdeTiIQDk4E2wFZgmYjMU9W1fs36APtUtY6IdAdGAzcDx4ChwIXuT1a3qGrWNUxzOpcpBnKb0ykiNYXBX8/0ladf2pEd5W0Zz4KQVCIKhg+HO+90Kl54Ae65B6pW9TYw44lAehSNgQRV3aiqycAsoHOWNp2BGe72XKC1iIiqHlbVb3ASRqCyPdcJHG9CVPdVn1Nr3w4AEqPK8HJTW8azQN1+O9Sv72wfOmRrVhRjgSSK6sCffuWtbl22bVQ1BUgEogM495vubaehfsngZM9lQljp5KM84Pdy3eTLb+JAybIeRlQMhIfDyJEZ5VdecR5um2LHy4fZt6jqRUAL9+e2EzlYRPqKSLyIxO+yeWlCXp9lH1Dl8H4Atper7NxDNwXv+uvhiiuc7ePHbWqPYiqQRLENqOlXruHWZdtGRCKACsCe3E6qqtvc3weBt3FucQV8LlWdqqqxqhpbpYrdpw5llY4k8n8/vucrj7viFuceuil4IjB6dEb5nXdgxQrv4jGeCCRRLAPqikgtEYkEugPzsrSZB/Ryt7sBi1VzfvdfRCJEpLK7XQK4HvjlZM5lQt/9382ibPJRANZXPov/XtjK44iKmWbNoEuXjHLcPwY+mhCXZ6JwnxP0BxYA64A5qrpGRIaLSCe32TQgWkQSgIcA398kEdkEjAV6i8hWEakHRAELRGQVsBKnF/FaXucyxU/N/X9xy0+f+srPX9WLtPQpJkzhee45CHM/LhYudH5MsRHQUqiq+gnwSZa6YX7bx4Bsh6CoakwOp700h/Y5nssUPwO/ecu3ct2PNeqxqHbjPI4wBeKCC5yhsq+/7pQfewyuucbeiC8m7M1sE7Tq7N5ClzVf+Mqjr+ptH0xeeuopKOmuQx4fD/Oy3oE2oSqgHoUxXnjwm7d9E/8tPieW5TXqeRxR8ZLdy4+PX9SOu5d94BSGDoWOHTNuSZmQZf+HTVCq9/dGrl//ja88tsWtHkZj0r3SpBuUcScMXL0a3n3X24BMobBEYYLSwG/+49v+7NzLbR3sILGnTEVn2dR0Tz4JKSneBWQKhSUKE3Qabl9Pm4QfAUhDGHfFLR5HZDIZPBjKl3e216+Ht97yNh5T4CxRmKDz0NcZvYmPLriS9VVivAvG/FOlSjBoUEb56adtcaMQZ4nCBJUmW1Zz5aafAGdRIptGPEg9+KCTMAD++APefNPbeEyBslFPptDlOJW4KrP9ehP/vbAVf1TKOv+kCQrly8Ojjzo/ACNGOLPNpg+fNSHFehQmaLTY9BNNtq4BIDksgonNe3gckclVv35w+unO9tatMHWqt/GYAmOJwgQHVQb59SZmX3wtWyuc7mFAJk9lyjhvaKd77jk4csS7eEyBsVtPJii0/v1HGu74DYCk8BK85L9mswkq/rcOo1JqsqRcZc48uBv+/hsmT4aHH/YwOlMQrEdhPCealqk38Z9G1/F3ucoeRmQClRQRyUvN/FYqHj0aDh70LiBTICxRGM9d+9sP1Nv5BwBHSkQxpWk3jyMyJ+Ldi65hS/ptwj17YNIkbwMy+c4ShfGUaBoP+i1xOrNRB3aXOc3DiMyJOh5egkn+vYoXX4QDB7wLyOQ7SxTGU+3Wf8cFuzYBcLhESaY26eptQOakvF+/FZxzjlPYtw8mTvQ2IJOvLFEYz4im8cC37/jKMy+5nr2lK3gYkTlZKeERmdfTHjMGEhO9C8jkq4AShYi0E5H1IpIgIv9YcU5EokRktrt/qYjEuPXRIrJERA6JyEt+7UuLyHwR+VVE1ojIKL99vUVkl4isdH/uOvXLNMHoul+/5fzdmwE4FFmKqY3/5XFE5pTceivUcSdv3L8fJkzwNh6Tb/IcHisi4cBkoA2wFVgmIvNUda1fsz7APlWtIyLdgdHAzcAxYChwofvj70VVXeKuw71IRNqravqal7NVtf8pXZnxXI5vYANhaamZehMzLrmefdabKNoi3F5FL3fJ+7FjnZlmK1b0Ni5zygLpUTQGElR1o6omA7OAzlnadAZmuNtzgdYiIqp6WFW/wUkYPqp6RFWXuNvJwAqgxilchyliOvz6Defu2QLAwchSvGa9idDQsyfUretsJybC+PHexmPyRSCJojrwp195q1uXbRtVTQESgehAAhCRikBHYJFfdVcRWSUic0WkZiDnMUVH1t7E9Es7sb9UeQ8jMvkmIgKGDcsojxvnPNw2RZqnD7NFJAJ4B5ioqhvd6o+AGFVtACwko6eS9di+IhIvIvG7du0qnIBNvui47ivq7N0KwIHI0rx+WRePIzL5qkcPOO88Z/vAAecWlCnSAkkU2wD/b/U13Lps27gf/hWAPQGceyqwQVV9/VNV3aOqSW7xdeDS7A5U1amqGquqsVWqVAngjzLBIDwtlQHfzfKV34ztRGKpch5GZPJdeHjmXsWECbB3r3fxmFMWyFxPy4C6IlILJyF0B7IuEjAP6AV8D3QDFquq5nZSEXkGJ6HclaW+mqrucIudgHUBxGiKiE5rv6T2Xud7xoGoMkyz3kTI8B+8EJZWlgXRNam7509nSo8xY+DZZz2MzpyKPHsU7jOH/sACnA/tOaq6RkSGi0gnt9k0IFpEEoCHAN8QWhHZBIwFeovIVhGpJyI1gMeBesCKLMNgB7hDZn8GBgC98+NCjffC01K53683MS22MwdKlvUwIlNQ0sLCmeA/TfzEibB7t3cBmVMS0OyxqvoJ8EmWumF+28eAG3M4NiaH00oO7YcAQwKJyxQtXdZ8wTn7tgOQGFWGN2M75XGEKcrmn38F9383i/N2b4FDh5xexciRXodlToK9mW0KRdbexOuXdbHeRIhTCWNCc7+71JMmgQ08KZIsUZhCccMvi4nZ7zx62l+yLG/GZn0Vx4SiT89rxroqMU7h8GFnwkBT5FiiMAUuIjUlU2/itcv+xaGo0h5GZAqL06vwe1bx0kuwc6d3AZmTYonCFLgbflnMWYl/A7CvZDlmXNrR44hMYVpw7uXQoIFTOHIEXnjB24DMCbNEYQpUidTj3P/9bF/5tcbWmyhuVMLgqacyKiZPdpZNNUWGJQpToLquXkRNtzext1R5ZlxyvccRGU906QINGzrbR4/C8897G485IZYoTIEpkXqc/n69iamNb+Cw9SaKJ5HMvYopU+CvvzwLx5wYSxSmwNy4+n/UOOAMh9xTqjwzL+ngcUTGU506wSWXONtHj8Lo0d7GYwJmicIUjKQk+n+X0Zt4tUlXjkSW8jAg47nsehXbt3sWjgmcJQpTMKZN48yDzpQNu0pX5N+NrDdhgOuvh9hYZzspCUaNyr29CQqWKEz+S0qC557zFV9p0pWjkSU9DMgEDRF4+umM8tSpsC3rZNQm2AQ015Mxucm65Gmv5R/xtPuPf1eZirzVqL0XYZlg1b49NGkCS5dmfKmYPNnrqEwurEdh8lXJ48cyjXSa0uRGjpWw3oTxk7VX8dprsHmzd/GYPFmiMPnqthWfUOXwfgB2lI223oTJ3rXXQvPmzvbx4zBihLfxmFzZrSeTb8okHeHepXN95Zea3UxSRKSHEZlgkfX2JEDTszow69tvncL06fDoo1C3buEGZgJiPQqTb3ov/4hKRw8A8GeF05nToI3HEZlg9sNZDaBVK6eQmpr5dpQJKtajMAHL7lthuvLHDtH3x/d85YnNunM8vERhhGWKshEjYPFiZ/vtt+Gxx6BePW9jMv8QUI9CRNqJyHoRSRCRuGz2R4nIbHf/UhGJceujRWSJiBwSkZeyHHOpiKx2j5koIuLWVxKRhSKywf192qlfpilofZZ9QIWkwwBsPO1M3ruwlccRmSKhWTNnFBSAKjz5pLfxmGzlmShEJByYDLTHWeO6h4hkTfl9gH2qWgcYB6S/m38MGAoMzubUU4C7gbruTzu3Pg5YpKp1gUX4rb9tgtNpRxLpE/+hrzz+ip6khoV7GJEpUvwfZM+dCytXeheLyVYgPYrGQIKqblTVZGAWkHV5ss7ADHd7LtBaRERVD6vqNzgJw0dEqgHlVfUHVVVgJtAlm3PN8Ks3Qer/fnyPsslHAfgt+iw+Pr+FxxGZIuXSS53ZZdMNG+ZdLCZbgSSK6sCffuWtbl22bVQ1BUgEovM459Ycznm6qu5wt/8CTg8gRuORKof20Wv5x77y2Ba3kGa9CROgmLj5xMTNp22la0hDnMqPPnJexjNBI6hHPbm9Dc1un4j0FZF4EYnfZQu2e+beH96lVEoSAGuqnuOsZmbMCVpfJYb551+RUWG9iqASSKLYBtT0K9dw67JtIyIRQAVgTx7nrJHDOf92b02l36LKdoFdVZ2qqrGqGlulSpUALsPktzMO7OaWlZ/6ymNa3OqsZmbMSRh/RU9S0//+fP45fPWVtwEZn0D+VS8D6opILRGJBLoD87K0mQf0cre7AYvd3kC23FtLB0SkqTva6XYg/Wmo/7l6+dWbIDPgu3eISj0OwE/VzmNx7cs8jsgUZb9H1+SD+ldnVDz2mDMSynguz0ThPnPoDywA1gFzVHWNiAwXkU5us2lAtIgkAA/hN1JJRDYBY4HeIrLVb8TUfcDrQALwO5D+1XQU0EZENgDXuGUTZM7Zs5WbVi30lce0uNWZw8eYUzC+eU+Opz/j+vZb+OQTbwMyQIAv3KnqJ8AnWeqG+W0fA27M4diYHOrjgQuzqd8DtA4kLuOdQV//mwhNA+Cbsy/mm1qNPI7IhII/K57B2w3b0WuF+3LnkCHQrh2E2wAJL9kNZXPCGuz4jQ7rv/WVn7+qVy6tjTkxk5p1hzJlnMLq1fDOO94GZCxRmBP3yJczfNvzz2vOqmrnehiNCTW7y5wGAwdmVAwdCsnJ3gVkLFGYE9N800qu2PwzACkSxpgWt3kckQlJgwdDtPsq1qZN8OqrnoZT3FmiMAETTePRL6f7ynMatGFjdI2cDzDmJMWM/IYRDTLe1t4dN5T6A9/1MKLizRKFCVj79d/R4K8EAI5FRDKheQ+PIzKh7D+XdGBbOecdqcpHErkz3kbKe8UShQnM8eMM/mqmrzj90o78Xa6yhwGZUJcUEcn4K3r6yn1/fA9sFgZPWKIwgXnzTc7Ztx2AA1FlmNKkm8cBmeLgvxe24rfoswAol3wURo70OKLiyRKFyduRI5lWH3ulSVcSS5XzMCBTXKSFhfPilX4DJiZPhi1bvAuomLJEYfI2aRJsd3oTf5etxBuxnfI4wJj883ndpqw48zynkJxsEwZ6wBKFyd3OnfDcc77ixGbdOVaipIcBmWJHhNFX9c4oz5wJP/3kWTjFka2Zbf7Bf23sEZ+/zG0HDgCQUKkGsxtc61VYphhbetZFLKzTmDYJPzoTBQ4aBIsW2fxihcR6FCZHdXdtpufKz3zl51reSUq4fbcw3hh59Z0Q4f79W7LEWeDIFAr7V29y9NgXbxDuN/GfTSNuvLQxugbTG7Sj9wpnRcWNve+lbZ80joeXAGDTqA5ehhfSrEdhstXijxW03LgcgDSEZ1v1sW6+8dyE5j04EOVMGHjOvu3c8tOneRxh8oMlCvMPYWmpPL54mq88p0Eb1lU9x8OIjHHsK12Bic1u9pUf/PZtKhw96GFExYMlCvMPN61ayPm7NwNwuERJZ1EiY4LEzEs6srniGQBUPHaI+7+b5XFEoc8Shcns4EEGff0fX3FK027sKlvJw4CMySw5ogQjr77DV759xXxi9m7zMKLQF1CiEJF2IrJeRBJEJC6b/VEiMtvdv1REYvz2DXHr14tIW7fuPBFZ6fdzQEQedPc9JSLb/PZdlz+XagIyahRVjuwHYHu5yrx+WZc8DjCm8H12bjN+rOGsqhyZlkKc36zGJv/lmShEJByYDLQH6gE9/Na9TtcH2KeqdYBxwGj32HpAd6A+0A54WUTCVXW9qjZU1YbApcAR4H2/841L3+8uw2oKw+bNMGaMr/j8Vb3s5ToTnER4ptVdvmK7376HL7/0MKDQFkiPojGQoKobVTUZmAV0ztKmM5C+7NlcoLWIiFs/S1WTVPUPIME9n7/WwO+quvlkL8LkkyFDICkJgJ/PqMuH9a7yOCBjcraq2rm8V79lRsVDD0FqqncBhbBAEkV14E+/8la3Lts2qpoCJALRAR7bHci6KG5/EVklIm+IyGnZBSUifUUkXkTid9nUw6fuyy8zrU38TKs+qNgjLBPcXrjydo5FRDqFFSvg9de9DShEefrCnYhEAp2AIX7VU4ARgLq/xwB3Zj1WVacCUwFiY2O1wIMNMf7TdESkpjB/+gDcadf4+PwWLKt5oTeBGXMCdpSvwpQm3Rj47dsA7H9wMC3XlGNf6Qq+NvYi3qkL5CvjNqCmX7mGW5dtGxGJACoAewI4tj2wQlXySdRNAAAWdklEQVT/Tq9Q1b9VNVVV04DX+OetKpPPei3/iPN2O1M3Hy5RkhGt+ngckTGBe6VJ10zDZR/5ckYeR5gTFUiiWAbUFZFabg+gOzAvS5t5QC93uxuwWFXVre/ujoqqBdQFfvQ7rgdZbjuJSDW/4r+AXwK9GHPiqh7c4/s2BjC+eU9buc4UKUklonjqmv/zlXus+pxG2371MKLQk2eicJ859AcWAOuAOaq6RkSGi0j6wgTTgGgRSQAeAuLcY9cAc4C1wGdAP1VNBRCRMkAb4L0sf+TzIrJaRFYBLYGBp3iNJhePL3mDsslHAfgt+izetLUmTBG0pPZlLKzTxFcevnAKYWn2YDu/BPSMwh2i+kmWumF+28eAG3M49lng2WzqD+M88M5af1vWOlMwLt+8is7rMoYUPtnmHpsd1hRZT7e+mxabfqJkSjIX/f07PX9ewH8a2WtY+cE+FUKc/0NrfxGpKQxfOMVX/vCCq/j+7AaFFZYx+W5rxTOY3PRGBn3zFgAPfzmDT85r7nFUocHGPxZTd8TPo+4eZ+TyochSPNvyHwPLjClypjbpyqaKzmPOCkmHefSL6d4GFCIsURRDZxzYzYN+D7DHNe/JznL/uAtoTJGTFBGZ6cH2zasXwvffexhRaLBEUQw9sWQaZY4fA2B95bOYcWlHjyMyJv98UTuWBXWbZlTcdx+kpHgXUAiwRFHMtEr4ket//dpXHtbmXnuAbULOiNZ3czQiyimsXAnjx3sbUBFnnxDFSLmkwzy7YLKv/N/6LVl61kUeRmRMwdha4XQmNbuZR76aCcCxIY/T7rfybKqUMYOQvbEdOOtRFCNDlrxBtUN7ANhVuiIjWt/tcUTGFJypjW9gjbsyY8mUZEZ/OhFx14A3J8YSRTHRbNNKev68wFce1uYe9pcq72FExhSslPAIHrnuAVLcyS2bbF1ja2yfJEsUxUDp5KOM+mySr/zpuc349PwrPIzImMKx5vTaTGma8S5w3JfTqZ6408OIiiZLFMXAw1/N5KxEZ97F/SXLMqzNvR5HZEzhmdSsOxuinblJyyYfZeRnk0BtwukTYYkixF26dS29ln/sKw9vfTe7yma7xIcxISk5ogSPth9AGgLAlZt+4sbV//M4qqLFEkUoO3aM5z+dSBjOt6cval3Ke/VbeRyUMYVvRfULeMNvwssnFr8O27d7GFHRYsNjQ0R2czo98uV07tu7FYCDkaV4rF0/ECns0IwJCi9eeRttEpZy9v6/qJB02HkR7/337d9EAKxHEaIabl9P36UZM7iPuvoOtpev6mFExnjrWImSxLUbkFHx4Yfw1lveBVSEWKIIQWWTjjDhoxeIcMeM/1DzQt5u2M7jqIzx3vdnN+Dti/3+Ldx3H2zc6F1ARYQlihA0fOEUzt7/FwAHIkszuMNAVOx/tTEAz7a80zfDLAcPQs+ecPy4t0EFuYA+PUSknYisF5EEEYnLZn+UiMx29y8VkRi/fUPc+vUi0tavfpO7kt1KEYn3q68kIgtFZIP724bonIAua5Zww5olvvLjbfuxtcLpHkZkTHA5HFWaAZ0e5nhYuFOxdCmTrrqNmLj5Oa7fUtzlmShEJByYDLQH6gE9RKRelmZ9gH2qWgcYB4x2j62Hs8Z2faAd8LJ7vnQtVbWhqsb61cUBi1S1LrDILZsA1Nz/FyM+f9lXfvfCa/io3lUeRmRMcFpV7VzGtMhYTLPf93NoumWVhxEFt0B6FI2BBFXdqKrJwCygc5Y2nYEZ7vZcoLWIiFs/S1WTVPUPIME9X278zzUD6BJAjMVeRGoKE+e9QDl3/es/TqvGU9f09TgqY4LXq01u4Ft3VccwlHEfjaHC0YMeRxWcAkkU1YE//cpb3bps26hqCpCIsx52bscq8LmILBcR/0+001V1h7v9F2D3TQLw4Ldv02jHegCOh4XzQMeHORxV2uOojAleKmEM7DCIve6cZ9UO7WH0ZxPtre1sePmE8wpVvQTnllY/EbkyawNVVSDb/2si0ldE4kUkfteuXQUcapD74gvu+/5dX/HFK29jVbVzPQzImKJhZ7loHmn/gK/c7rfv4bXXPIwoOAWSKLYBNf3KNdy6bNuISARQAdiT27Gqmv57J/A+Gbek/haRau65qgHZzuClqlNVNVZVY6tUqRLAZYSoPXvg1lt9b19/c/bFTG18g8dBGVN0/K9uE2Y28lub4sEHYe1a7wIKQoG8mb0MqCsitXA+5LsDPbO0mQf0Ar4HugGLVVVFZB7wtoiMBc4E6gI/ikgZIExVD7rb1wLDs5xrlPv7w1O5wFDjPyojPC2VN999iiu3OXl7b6nyPNThIRsKa8wJerblnTT5czXn7d4CR4/ye4u2dLl9DAejyvjaFOeFjvL8RHGfOfQHFgDrgDmqukZEhotI+uQp04BoEUkAHsIdqaSqa4A5wFrgM6CfqqbiPHf4RkR+Bn4E5qvqZ+65RgFtRGQDcI1bNtl4+MsZXLnpp4zydQ+ws1y0hxEZUzQllYhiQKdHfMun1t67lbEfj7GFjlyiIfDgJjY2VuPj4/NuGALSexSd1n7JxI9e8NVPaNadcS1u9SosY0JCp7VfMPGjF33lCc16MK7FLUBo9ihEZHmW1xOyZfcoiqD6f//O6E8n+sr/q30Z46/IejfQGHOi5tW7mlf9nvE98N07tP3tOw8jCg6WKIqYSkcSefW9ZyiVkgTA75VqMLDjYHsuYUw+ef6qXnwV08hXHjN/HHV3bfYwIu/Zp0tRcvw4kz8cRY0DznDgA5GlufuGJzI9cDPGnJrUsHDu7/QImyueATir4k19/xnYt8/jyLxjiaIoGTyYy7esBiAN4cGOg9kYXcPjoIwJPYmlytH3hic4XKIkALX27XAmD0xN9Tgyb9jCRUEq6+RkN/38Oc9/lvFcYkyLW1lcJ6/ZUIwxJ2t9lRgGdRjIKx+MdCo++wzi4uCFF3I/MARZj6IIuGbDUp5b8JKv/Mm5zZh8+U0eRmRM8fDZec2ZdPnNGRUvvghjx3oXkEcsUQS5y/78hZfmjfYtQrS2ai0GdxhoyzcaU0jGtriFz+s2zagYNAj+/W/vAvKAJYogdt6uTUz77whKpiQDsLniGfS6cThHIkt5HJkxxYdKGPd3fJgfa2SsrnC89x30vvHpYrOGhSWKIFUj8W9mzhlG+aTDAOwqU5HbbhrBrrK2jpMxhS2pRBR3dR3GuioxAJRIS2XKByO5ZNs6bwMrJJYogtHOncycPZTTD+0FnGGwvW98mi2nVfM4MGOKrwMly9Lrxqf5010xslRKEm/MfZo6u7d4HFnBs0QRbA4ehOuu45x92wFICi9B365PsOb02h4HZozZWS6a224azu7SFQCoeOwQM+cMgy2hnSwsUQSTgwehY0dYvhyAVAljQMeH+eGsBh4HZoxJt6lSdXrf+DSH3GeFZx7cDW3awNatHkdWcOw9Cg/5PwQrf+wQ0999kku2r/fVPXHtfSw4r5kXoRljcvHLGXXo+6/HeXPuU0SlpsBvv0GLFrBoEZxzjtfh5TvrUQSBSkcSeeedxzIlieeuvoN3GrbzMCpjTG6+i2nIgE6PkBzmft/etAmuvBJ+/dXTuAqCJQqPVT24h1lvD6H+zo2+umHX/B9Tm3T1MCpjTCAWnNuMvjc8DlHOOhZs2+Yki1WrvA0sn1mi8FD1xJ3MeTuOc/c4D8JSJYyH2z/AzEs7ehyZMSZQX9S+DD75BMq4k3Pu2gVXXw3LlnkaV36yROGVDRuY89ajxOzfAcDxsHAe6DiYdxu08TgwY8yJivn8KDd0eYoD6TM579vHwSuugq+/9jawfBJQohCRdiKyXkQSRCQum/1RIjLb3b9URGL89g1x69eLSFu3rqaILBGRtSKyRkQe8Gv/lIhsE5GV7s91p36ZQebbb+GKK6h+0JkuPCk8gnu7PMbHF1zpcWDGmJO1osYF9Oj+LHtLlQegXPJRuPZamD3b48hOXZ6JQkTCgclAe6Ae0ENE6mVp1gfYp6p1gHHAaPfYekB3oD7QDnjZPV8KMEhV6wFNgX5ZzjlOVRu6P5+c0hUGmzffhJYtYedOAI5GRNGn65P8r24TjwMzxpyqNWfU4eYeI9lZxp1B4dgx6N4dhg6FtKK7/nYgPYrGQIKqblTVZGAW0DlLm87ADHd7LtBaRMStn6WqSar6B5AANFbVHaq6AkBVDwLrgOqnfjlBLCUFHnoI7rwTjh8HYE+p8tx68wi+qdUoj4ONMUXFhipnc+Mto/m9kt9H2jPPQLducOiQd4GdgkASRXXgT7/yVv75oe5ro6opQCIQHcix7m2qRsBSv+r+IrJKRN4QkWwnNxKRviISLyLxu3btCuAyPLRvH3ToAOPGZdQ1aEDnXuNYXiNr58wYU9RtPu1M/nXbGGjbNqPy/fehWTNnGG0R4+nDbBEpC/wXeFBVD7jVU4DaQENgBzAmu2NVdaqqxqpqbJUqVQol3pOyfj2/17kIPv/cV/XZuZdTr/VQtrpzxhhjQs+BkmWp3eA+XrusS0bl6tXsqXcxfPWVd4GdhEASxTagpl+5hluXbRsRiQAqAHtyO1ZESuAkibdU9b30Bqr6t6qmqmoa8BrOra+iRxXeegsaN6b23oz/XBOadefeLkNsqnBjioHUsHCebXUXD7d/gKRw58W86KMHoHVrGD26yCytGkiiWAbUFZFaIhKJ83B6XpY284Be7nY3YLGqqlvf3R0VVQuoC/zoPr+YBqxT1UzLRYmI/xSp/wJ+OdGL8ty+fdCjB9x6KxxwOkpHI6Lo1+lRxrW4FRUblWxMcfJugzb06D6SXaUrOhUpKc6yqi1bFolbUXl+YrnPHPoDC3AeOs9R1TUiMlxEOrnNpgHRIpIAPATEuceuAeYAa4HPgH6qmgo0B24DWmUzDPZ5EVktIquAlsDA/LrYQrFoEVx0UaYhcZsrnkG3W0Yz/4IWHgZmjPHSihoX0LnXWFZWOzej8uuvoUEDmDnTuQsRpESDOLhAxcbGanx8vLdBHDvG6y26c1f8h5mqZzW4lhGt7uJwVGmPAjPGBJOI1BQSyvzkjITyv/XUrRu88gpERxdaLCKyXFVj82pn90Dyw3ffwWWXZUoSe0qVp++/Hieu/QBLEsYYn5TwCHjqKefF2zp1MnbMnevcjXjvvaDrXViiOBXbt8Ntt0Hz5vBLxqOUxefE0u7OyXx+7uUeBmeMCVYxcfOJeX839TqO5O2L/WaJ3rEDunZ13uheFzzLrFqiOAl1H/6AkS3v4FCtOvCf//jqj0ZE8fi193FntydtbWtjTJ6ORJbisXb96dN1aMaDboD//c95dvHQQ5CY6F2ALntGcaI+/ZTfb70705BXgPnnNee5ln3YVqFq4cRhjAkpFY4e5OfUb2DKlMzTfVStCqNGQa9eEJa/3+3tGUV+UnVGM119NVx3XaYksb7yWfTo/iz9ugyxJGGMOWmJpcoRU7Y97XpN4IeaF2bs2LnTmfqnYUNnNKUH715YjyIXMY9+zNUbl3P/d7O4dHvmVasORJVh7BW38J9G1zkPp4wxJr+ocv2vX/PYkjecNbn9nX8+PPaY865WxKl99gTao7BEkZ20NJg3j1X3PkyDvxIy7ToeFs67F13DmBa3sadMxRxOYIwxp65U8jHu/eFdBqz6GA4fzryzdm0YMsQZUBMZeVLnt0RxMnbtgunTYepUSMicIJLCI3j3oja80rSbzdFkjClUmwY3gfHjYdIk32wPPmedBS+9BB1PfGXMQBNFsb9nEvPox1y+ZTU9V35K29++JzItJdP+YxGRvHNxW15t3JW/ylf2KEpjTHEW8+JS4HLK33ERvZd/xJ3xH1LxmDtl+ZYtGcuwFpDi26Nwew8bR47nnH3b/7H7QGRpZl3cltca32BDXY0xQaVs0hFu/ekT7lr2PpUb1odvvgGREz6P9Sjy8sYbEBfHOVmqV5x5Hu9c3I6Pz2/B0ciSnoRmjDG5ORRVmleadmPGJdezrn+jk0oSJ6L4Joo77nCWJzx+nAORpXn/wpbMurgt66pmTR3GGBOcjkaWhHMK/jOr+CaKqlVhyBAeXpZovQdjjMlF8X7h7umnebdBG0sSxhiTi+KdKIwxxuTJEoUxxphcBZQoRKSdiKwXkQQRictmf5SIzHb3LxWRGL99Q9z69SLSNq9zukuuLnXrZ7vLrxpjjPFInolCRMKByUB7oB7QQ0TqZWnWB9inqnWAccBo99h6OGts1wfaAS+LSHge5xwNjHPPtc89tzHGGI8E0qNoDCSo6kZVTQZmAZ2ztOkMzHC35wKtRUTc+lmqmqSqfwAJ7vmyPad7TCv3HLjn7HLyl2eMMeZUBZIoqgN/+pW3unXZtlHVFCARiM7l2Jzqo4H97jly+rOMMcYUoiL7HoWI9AX6usVDIrLey3hOUmVgd56tQk9xvW4ovtdu111AZPQpHX52II0CSRTbgJp+5RpuXXZttopIBFAB2JPHsdnV7wEqikiE26vI7s8CQFWnAlMDiD9oiUh8IPOshJriet1QfK/drrtoC+TW0zKgrjsaKRLn4fS8LG3mAb3c7W7AYnVmG5wHdHdHRdUC6gI/5nRO95gl7jlwz/nhyV+eMcaYU5Vnj0JVU0SkP7AACAfeUNU1IjIciFfVecA04N8ikgDsxfngx203B1gLpAD9VDUVILtzun/ko8AsEXkG+Mk9tzHGGI+ExDTjRZWI9HVvoRUrxfW6ofheu1130WaJwhhjTK5sCg9jjDG5skThIREZJCIqIpXdsojIRHf6klUiconXMeYnEXlBRH51r+19Eanoty/bqV5CRV7T4IQKEakpIktEZK2IrBGRB9z6SiKyUEQ2uL9DctlId+aJn0TkY7ccElMSWaLwiIjUBK4FtvhVt8cZGVYX5x2RKR6EVpAWAheqagPgN2AI5DzVi2dR5rMAp8EJFSnAIFWtBzQF+rnXGgcsUtW6wCK3HIoeANb5lUNiSiJLFN4ZBzwC+D8k6gzMVMcPOO+UVPMkugKgqp/7vXX/A857MpDzVC+hIpBpcEKCqu5Q1RXu9kGcD83qZJ7mJySn5hGRGkAH4HW3HDJTElmi8ICIdAa2qerPWXYFMl1KqLgT+NTdDvXrDvXry5Y7i3QjYClwuqrucHf9BZzuUVgFaTzOl780txwyUxIV2Sk8gp2I/A84I5tdjwOP4dx2Cjm5Xbeqfui2eRznFsVbhRmbKTwiUhb4L/Cgqh5wvlw7VFVFJKSGW4rI9cBOVV0uIld7HU9+s0RRQFT1muzqReQioBbws/uPpwawQkQaE9h0KUEtp+tOJyK9geuB1poxNrvIX3ceQv36MhGREjhJ4i1Vfc+t/ltEqqnqDvd26k7vIiwQzYFOInIdUBIoD0wgwCmJgp3deipkqrpaVauqaoyqxuB0Ry9R1b9wpjy53R391BRI9OuuF3ki0g6na95JVY/47cppqpdQEcg0OCHBvS8/DVinqmP9dvlP8xNyU/Oo6hBVreH+m+6OM43RLYTIlETWowgunwDX4TzMPQLc4W04+e4lIApY6PamflDVe3Kb6iUU5DQNjsdhFZTmwG3AahFZ6dY9BowC5ohIH2AzcJNH8RW2kJiSyN7MNsYYkyu79WSMMSZXliiMMcbkyhKFMcaYXFmiMMYYkytLFMYYY3JlicIYY0yuLFEYY4zJlSUKY4wxufp/FAnOSE0Vm34AAAAASUVORK5CYII=\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 = 20.\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",
    "#    weight = np.exp(-0.5*(xtrial-xold)/sigma2)\n",
    "#    if(xtrial < 0):\n",
    "#        weight = 0\n",
    "\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 = 5\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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}