tatsy/hmc.ipynb

## hmc.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "n_mutate = 10000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inv_dist_fun(x):\n",
    "    r = 0.3\n",
    "    mu1 = 1.0\n",
    "    mu2 = -2.0\n",
    "    x1 = x - mu1\n",
    "    x2 = x - mu2\n",
    "    e1 = np.exp(-x1 * x1) / np.sqrt(np.pi)\n",
    "    e2 = np.exp(-x2 * x2) / np.sqrt(np.pi)\n",
    "    return r * e1 + (1.0 - r) * e2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class HMC(object):\n",
    "    def __init__(self, x, inv_dist_fun, sigma, seed=0):\n",
    "        self.rng = np.random.RandomState(seed)\n",
    "        self.x = x\n",
    "        self.inv_dist_fun = inv_dist_fun\n",
    "        self.sigma = sigma\n",
    "        self.L = 50\n",
    "        self.delta = 0.03\n",
    "\n",
    "    def sample(self):\n",
    "        x0 = self.x\n",
    "        u0 = self.rng.normal(0.0, 1.0)\n",
    "        eps = 1.0e-12\n",
    "\n",
    "        # Leap-frog\n",
    "        x1 = x0\n",
    "        u1 = u0\n",
    "        for l in range(self.L):\n",
    "            x_temp = x1 + 0.5 * self.delta * u1\n",
    "            x_plus = np.log(self.inv_dist_fun(x_temp + eps))\n",
    "            x_minus = np.log(self.inv_dist_fun(x_temp - eps))\n",
    "            dx = (x_plus - x_minus) / (2.0 * eps)\n",
    "            u1 = u1 + self.delta * dx\n",
    "            x1 = x_temp + 0.5 * self.delta * u1\n",
    "\n",
    "        exp_H0_minus = self.inv_dist_fun(x0) * np.exp(-0.5 * u0 * u0)\n",
    "        exp_H1_minus = self.inv_dist_fun(x1) * np.exp(-0.5 * u1 * u1)\n",
    "        a = min(1.0, exp_H1_minus / (exp_H0_minus + 1.0e-12))\n",
    "        if self.rng.random() < a:\n",
    "            self.x = x1\n",
    "        else:\n",
    "            self.x = x0\n",
    "\n",
    "        return self.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler = HMC(0.0, inv_dist_fun, 1.0)\n",
    "\n",
    "# Burn-in\n",
    "n_burnin = max(n_mutate // 100, 100)\n",
    "for i in range(n_burnin):\n",
    "    sampler.sample()\n",
    "\n",
    "# Sampling\n",
    "samples = []\n",
    "for i in range(n_mutate):\n",
    "    samples.append(sampler.sample())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAulUlEQVR4nO3deXxU1d3H8c+ZZLKvJAQSkpCAYZMAQlgURC3iA4igFitYra21aFtqxWrFpa1V28e1WitKqbti8VGrIqCioCKCQNgJEAyQhJAAWcky2WbmPH8k0AABJsnM3Fl+79fLl5m5d+75TjS/nJw59xyltUYIIYT3MxkdQAghhHNIQRdCCB8hBV0IIXyEFHQhhPARUtCFEMJHSEEXQggf4VBBV0pNUkrlKqXylFLz2jl+qVLqmFJqa+s/f3R+VCGEEGcTeK4TlFIBwHxgIlAEbFRKLdFa7zrl1G+01lNdkFEIIYQDHOmhjwLytNb7tdZNwGJgumtjCSGE6Khz9tCBXsDBNo+LgNHtnHehUmobUAzcrbXOOfUEpdRsYDZAeHj4iAEDBnQ8sRBC+LFNmzaVaa27t3fMkYKu2nnu1PUCNgO9tda1SqkpwIdAxmkv0nohsBAgKytLZ2dnO9C8EEKI45RSBWc65siQSxGQ0uZxMi298BO01tVa69rWr5cDZqVUfCeyCiGE6CRHCvpGIEMpla6UCgJmAkvanqCU6qmUUq1fj2q9brmzwwohhDizcw65aK2tSqk5wGdAAPCK1jpHKXV76/EFwAzgl0opK1APzNSyjKMQQriVMqruyhi6EP6lubmZoqIiGhoajI7iFUJCQkhOTsZsNp/0vFJqk9Y6q73XOPKhqBBCdFlRURGRkZGkpaXROkIrzkBrTXl5OUVFRaSnpzv8Orn1XwjhFg0NDcTFxUkxd4BSiri4uA7/NSMFXQjhNlLMHdeZ75UUdCGE8BFS0IUQfqGqqooXXnjB5e18+OGH7Np16lJX7iEFXQjhFzpa0LXW2O32DrcjBV0IIVxs3rx57Nu3j2HDhjF37lwmTJjA8OHDyczM5KOPPgIgPz+fgQMH8qtf/Yrhw4dz8OBBHnnkEQYMGMDEiROZNWsWTz31FAD79u1j0qRJjBgxgosvvpg9e/awdu1alixZwj333MOwYcPYt2+fW9+jTFsUQhii4KafnPZc5ORJdLvhBuz19Rycfdtpx6OvuYaYa6/BWlnJoTt+e9Kx3m++cdb2HnvsMXbu3MnWrVuxWq1YLBaioqIoKytjzJgxTJs2DYDc3FxeffVVXnjhBbKzs3n//ffZsmULVquV4cOHM2LECABmz57NggULyMjIYP369fzqV79i1apVTJs2jalTpzJjxozOfms6TQq6EMLvaK25//77Wb16NSaTiUOHDnHkyBEAevfuzZgxYwBYs2YN06dPJzQ0FICrrroKgNraWtauXct111134pqNjY1ufhenk4IuhDDE2XrUptDQsx4PjI09Z4/8bBYtWkRpaSmbNm3CbDaTlpZ2Ys53eHj4ifPOdCe93W4nJiaGrVu3djqDK8gYuhDCL0RGRlJTUwPAsWPHSEhIwGw28+WXX1JQ0P6KtOPGjePjjz+moaGB2tpali1bBkBUVBTp6em8++67QEvh37Zt22ntuJsUdCGEX4iLi2Ps2LEMHjyYrVu3kp2dTVZWFosWLeJMm+2MHDmSadOmMXToUK699lqysrKIjo4GWnr5L7/8MkOHDuX8888/8cHqzJkzefLJJ7ngggvc/qGoLM4lhHCL3bt3M3DgQKNjdFhtbS0RERFYLBbGjx/PwoULGT58uFvabu97JotzCSFEJ82ePZtdu3bR0NDAzTff7LZi3hlS0IUQ4izefvttoyM4TMbQhRDCR0hBF0IIHyEFXQghfIQUdCGE8BHyoagQwhBp85Y59Xr5j13p1Ot1xaWXXspTTz1FVla7swtdRnroQgjhI6SgCyH8Ql1dHVdeeSVDhw5l8ODBvPPOOzz88MOMHDmSwYMHM3v27BNrt1x66aXMnTuX8ePHM3DgQDZu3Mi1115LRkYGDz74INCy1O6AAQO4+eabGTJkCDNmzMBisZzW7ooVK7jwwgsZPnw41113HbW1tUDLcr6DBg1iyJAh3H333U55j1LQhRB+4dNPPyUpKYlt27axc+dOJk2axJw5c9i4cSM7d+6kvr6epUuXnjg/KCiI1atXc/vttzN9+nTmz5/Pzp07ee211ygvLwdaltqdPXs227dvJyoq6rQNNMrKynj00Uf54osv2Lx5M1lZWfztb3+joqKCDz74gJycHLZv337il0RXSUEXQviFzMxMvvjiC+69916++eYboqOj+fLLLxk9ejSZmZmsWrWKnJycE+cfXx89MzOT888/n8TERIKDg+nTpw8HDx4EICUlhbFjxwJw4403smbNmpPa/O6779i1axdjx45l2LBhvP766xQUFBAVFUVISAi33nor//nPfwgLC3PKe5QPRYUQfqFfv35s2rSJ5cuXc99993HFFVcwf/58srOzSUlJ4aGHHjqxhC5AcHAwACaT6cTXxx9brVYAlFIntXHqY601EydO5N///vdpeTZs2MDKlStZvHgxzz//PKtWrerye5QeuhDCLxQXFxMWFsaNN97I3XffzebNmwGIj4+ntraW9957r8PXLCwsZN26dQD8+9//Zty4cScdHzNmDN9++y15eXkAWCwW9u7dS21tLceOHWPKlCk8++yzTltXXXroQghDuHua4Y4dO7jnnnswmUyYzWZefPFFPvzwQzIzM0lLS2PkyJEdvubAgQN5/fXXue2228jIyOCXv/zlSce7d+/Oa6+9xqxZs07saPToo48SGRnJ9OnTaWhoQGvNM88845T3KMvnCiHcwluXzz2T/Px8pk6dys6dO13WRkeXz5UhFyGE8BFS0IUQohPS0tJc2jvvDCnoQgi3MWqI1xt15nslBV0I4RYhISGUl5dLUXeA1pry8nJCQkI69DqZ5SKEcIvk5GSKioooLS01OopXCAkJITk5uUOvkYIuhHALs9lMenq60TF8mgy5CCGEj3CooCulJimlcpVSeUqpeWc5b6RSyqaUmuG8iEIIIRxxzoKulAoA5gOTgUHALKXUoDOc9zjwmbNDCiGEODdHxtBHAXla6/0ASqnFwHRg1ynn/QZ4H+j4/bPCI5xpBxlP2glGCHFmjgy59AIOtnlc1PrcCUqpXsA1wIKzXUgpNVspla2UypZPuoUQwrkc6aGrdp47dSLps8C9WmvbqctHnvQirRcCC6FlLRcHMwon68hejqHNDZT8+c/ETJ9O6LBhrgslhOgyR3roRUBKm8fJQPEp52QBi5VS+cAM4AWl1NXOCCjcr0/VIS4q3gFAYl051R8tIX/mLI4+/TTabjc4nRDiTBzpoW8EMpRS6cAhYCZwQ9sTtNYnJpcqpV4DlmqtP3ReTOEuibVl/HXtQmqCwljfcxD7Y3qR8c1qjjz+BOX/eglMASTMvdPomEKIdpyzoGutrUqpObTMXgkAXtFa5yilbm89ftZxc+E9Auw27tv4JqD545ifYzMFAGAKD6fnnx8CrSn/5z8JHTaUyMsuMzSrEOJ0Dt0pqrVeDiw/5bl2C7nW+qddjyWM8MO8r8g4doi/jLyJkoj4k44ppej5h5aNbIP79DEinhDiHOTWfwFAdGMNM3NXsjZxMGt6DW33HBUUROIjD7s5mRDCUXLrvwCgJiic+UOv5ZXzzz3nvKnoEMX3zsNaUeGGZEIIR0kPXQBgVyZWpra7q9Vp0xxTao6wYOUSzL2S6H7HHe6IJ4RwgBR0wf/kryeqqY53My6Ds9xHcNzByB5s7DGA/i+9wUUHU2kOMJ90XO4sFcIYMuTi50zazqzczxl+NNehYn7ch+eNJ6apjkuLtrgwnRCiI6Sg+7lRh3fRo76KpeljO/S6rfHnURiZwP8UbHBRMiFER8mQi5+7omAj5SFRrEs8v2MvVIoP+1xMSu1RTNqOXUnfQAijSUH3YxFNFrKO7GFJn7HYW28i6ohP0i90QSohRGdJt8qPRTZb2JzQj69Shnf6GkrbGVy2D6VljRchjCYF3Y+VhMfz0IU/Jy+mYxvRtjX+0DaeXPMi/SoPnvtkIYRLSUH3U2HNDcTXV3X5OpsS+mNVJi4q2dn1UEKILpGC7qfGFW/nzc8eJbnmaJeuUxsUxvb481oKupYl7oUwkhR0PzW6JIejoTEURXTv8rXWJg0mubaU1JojTkgmhOgsKeh+yGxrZnjpXtb3HNShm4nOZF3PlimPF5bkdPlaQojOk2mLfmhYaR4htma+69nBuednUBEazW8vuYP90UlOuZ4QonOkoPuh0Yd3YQkMZkd8X6ddc29sqtOuJYToHCnofmjRgIl8m5RJc4Dz/vOHNjcwa+8XbOneD5DFuYQwghR0P1QZEkVlSJRTr9kYYGZS/nqiGi1Ova4QwnHyoaifqf3mG6bu/xaTk+/stJsC2NI9gxFH96Bl+qIQhpCC7mcq33mHH+Z9hZ2uz2451aaE/sQ3VNP4/fdOv7YQ4tykoPsRbbNhWb+Brd0znDJd8VSbe/QHoG7Nt06/thDi3KSg+5GGnBzsNTUtBd0FykJj2B3bG93U5JLrCyHOTj4U9SN1674DYGv381zWxl2X/Ib822WWixBGkB66H2kuKSZ44ECOBUe6vC35YFQI95OC7kcSH3qI9HcWu7SNYGsj+6ZcScVrr7u0HSHE6aSg+xkVFOTS6zcGBoPdjmWD7DUqhLtJQfcTFW+8SeGtv0A3N7u8rbBRo7BkZ6NtNpe3JYT4LynofqJ2zTc0Hy5Bmc0ubyts1CjsNTU07Nnj8raEEP8lBd0PaJuN+i1bCRs+wi3thY0cCYBlvQy7COFOUtD9QGNeHvaaGsKy3FPQzT0S6PbzWwju388t7QkhWsg8dB+VNm/Zia+n7v+WXwMTPz/G0bXLzvwiJ+pxzz1uaUcI8V/SQ/cDVcERfJM0hKNhsW5rU2tN44ED2I4dc1ubQvg76aH7gTW9hrKm11C3tZc2bxm9q0tYsOppnho+k5WpWQDkPyZ3kArhStJD93HB1iZCrI1ub7cwsgd1gSEMrMh3e9tC+CuHCrpSapJSKlcplaeUmtfO8elKqe1Kqa1KqWyl1DjnRxWdMbZ4B+8t+wNJtaVubVcrE7u79WZQRYFb2xXCn52zoCulAoD5wGRgEDBLKTXolNNWAkO11sOAW4CXnJxTdNLAynwaA8wcDo9ze9u7u/Wmd/Vhwprr3d62EP7IkR76KCBPa71fa90ELAamtz1Ba12r/7saUzggKzN5iH6VB9kbk4JduX90bVdcOiY0AyoL3d62EP7IkZ/yXsDBNo+LWp87iVLqGqXUHmAZLb300yilZrcOyWSXlrp3CMAfmW3N9DlWTG5sqiHt745N5ZFRN7M3JsWQ9oXwN44U9Pa2tjmtB661/kBrPQC4GnikvQtprRdqrbO01lndu3fvUFDRcX2PFROo7eyNNaagNgYGszYpk9qgMEPaF8LfOFLQi4C2FSEZKD7TyVrr1UBfpVR8F7OJLioNjeGfg6eRE5dmWIbE2jKuyfva6ZtSCyFO50hB3whkKKXSlVJBwExgSdsTlFLnKdWySaVSajgQBJQ7O6zomPLQaD48b7xbNrQ4kwGVhcze+TG9qw8blkEIf3HOG4u01lal1BzgMyAAeEVrnaOUur31+ALgh8BPlFLNQD1wvZYtaww38vBu8mJ6URkSZViGXd16AzBI5qML4XIO3SmqtV4OLD/luQVtvn4ceNy50URXRDRZePi7l3l10GT+r98Ew3IcCetGVVA4/SoPnvtkIUSXyJ2iPqpfVUsB3WvQDJcTlCI3NpX+MnVRCJeTgu6jjveI98YkG5wE9samkFRXhr1ebjASwpVkcS4f1b+ykMKIBCzmUKOj8GHfi3kv4zL2hhqfRQhfJj10H6S1brlD1KD556eymENpCnD91ndC+Dvpofuou8f/GuVBE42uzltN6XP76X7Hb4yOIoTPkh66D1JKURIeT3GE59yN26/qIFXvvWd0DCF8mhR0H1T9ySdckb/e6Bgn2RObivXoUZqPHDE6ihA+Swq6D6p8+99MLvjO6BgnOT59sn77doOTCOG7pKD7GG2zUZ+TQ25sb6OjnGRfdBKYzTRIQRfCZaSg+5jGvH1oi4VcD5nhclxzgJnwUaOMjiGET5NZLj6mYUdLD9jwO0TbkfqybGQlhCtJD93HNBUVYYqOptiALeeEEMaSgu5jEu68k4yvvkQbsOXcuVjLyth/1VVUffCh0VGE8Eme91MvuszkobfYB3TrRvPhI9Rv3Wp0FCF8khR0H1K/YyeFt91G44EDRkdplzKZCM0cTP0OmekihCtIQfch9Zs3Uff1akzh4UZHOaOQzCE05u6VlReFcAGZ5eJD6rfvIDAxEXNCgtFR2pU2bxmjS5p4yGbjqjteJicuHYD8x640OJkQvkEKupdLm7fsxNevfPkd+2KSmNjmOU+TG5vKyuTh1AcEGR1FCJ8jQy4+IqqxjkRLObkeOP+8raqQSJ7KuoH9Mb2MjiKEz5GC7iMim+rYEdeH3d3SjI5yblrT3VJpdAohfI4UdB9xKDKB31/8qxPj0p5s2v41vLHiL0Q31hodRQifIgXdR5i03egIDtsX3TLc0k82jhbCqaSg+wKtef2zR7lhzwqjkzgkL6YXNhQDpKAL4VRS0H1AD0sF8Q3VVAVHGB3FIY2BwRRE9aRf5UGjowjhU6Sg+4D+rYXR02e4tJUbm0r/ykLwoH1PhfB2Mg/dB/SrLKTRFEh+VKLRURz2We9RbEnI8KqxfyE8nRR0H9C/6iB5McnYTAFGR3FYbrfe5OJZuyoJ4e2koPuAb5KGYgkMNjpGh51XVUSotdHoGEL4DCnoPmBJ33FGR+iUX+xYQpC9GbjD6ChC+AT5UNTLdbdUEtNQY3SMTsmNTaXvsWLsTU1GRxHCJ0hB93I35H7OP1c+6ZWzRXJjUzHbbTTu2WN0FCF8ghR0L9e/8iB7Y1NAKaOjdNjxaZb122TDCyGcQQq6F7PX1ZFafdir5p+3VRYaTXlIFPXbpaAL4QxS0L1Y/c4cAtBeW9BRigcu+gU9//Qno5MI4ROkoHux+u3bAO+6Q/RUBVGJBER47pZ5QngThwq6UmqSUipXKZWnlJrXzvEfK6W2t/6zVik11PlRxamip0zhryNvpDrYewtiZFMdpc89R/2OnUZHEcLrnXMeulIqAJgPTASKgI1KqSVa611tTjsAXKK1rlRKTQYWAqNdEVj8l7lXL77pNczoGF1iUybKXlwAgYGEZg42Oo4QXs2RHvooIE9rvV9r3QQsBqa3PUFrvVZrfXwLmu+AZOfGFKeyVlZS9d57RDd65xz04yzmUIL69KFh+w6jowjh9Rwp6L2AtuucFrU+dyY/Bz5p74BSarZSKlsplV1aWup4SnEaS3Y2JQ/+gcS6CqOjdFnokCHUb9+O9sK59EJ4EkcKensTnNv9yVNKXUZLQb+3veNa64Va6yytdVb37t0dTylO07B9O5jN7ItOMjpKl4UOycRWUUHzoWKjowjh1Rwp6EVASpvHycBpP3lKqSHAS8B0rXW5c+KJM6nftp2QAQNoDjAbHaXLQoYMwRQWRnORbHghRFc4UtA3AhlKqXSlVBAwE1jS9gSlVCrwH+AmrfVe58cUbWmbjYadOwnNzDQ6ilOEDBxIv40bCB8zxugoQni1c85y0VpblVJzgM+AAOAVrXWOUur21uMLgD8CccALquUWdKvWOst1sf1bU0EBdouF0KFDYJ3RabpOmeR2CCGcwaHlc7XWy4Hlpzy3oM3XtwK3OjeaOJPgPn3IWPstKigY1n1ldBynqP78c8oX/ovei97CFBRkdBwhvJJ0jbxUYLduPneHZcOOHTTszDE6hhBeSwq6Fzr8179S/Um7M0O9VtgFFwBQv2WzwUmE8F5S0L2M3WKh8q1FNH7/vdFRnCowPh5z71Qsm7cYHUUIryUF3cs05OSA3U7IkCFGR3G6sAuGU795s9xgJEQnyZ6iXqZ+W8sKi6FDfW/9s4hLxqObGrHXWXzu8wEh3EEKupexbN5CUFoagbGxRkdxuqjJk4maPNnoGEJ4LSno3kZrwsb41kKWafOWnfQ4rLkeizmU/MeuNCiRezUVHaJu7bc0FxYSMWHCiQ+IhegoKeheJuXFF3x6jPl3m/5N/8pCZl/e7nJAPqWpoIAjTzxJ7apVoDXKbCZsdMvdspbNmyl/+RUS7ppLcN++BicV3kI+FPVCygs3hHbUwcgEUmpLiWqsMzqKS9lq6zhw3Y+wrF9P/C9vp8/y5fTfvo2Ii8cB0HyoGMuGDRy45loq3lrk07/EhfNID92LHH74EZqLi0lZ8KLRUVxmV7c0AAZW5Buaw9UCIsJJfPQRQocOw9wj4bTj0VdNJfzCMZQ8+AeOPPoozQcLSbj3XlkmQZyVFHQvkTZvGf9c+QUlYXE8dMqYsy/ZG5tKswrg/IoDRkdxOq01R/73fwkdMpToqVcSdcUVZz0/MD6e5Bfmc/Txx6l4/Q2CBw4k5uqr3RNWeCX5de8lIpospNYcZVdcmtFRXKopwMy+mF4MKs83OorTlT0/n8o33qRxz26HX6NMJnrcdx/JC14keto0F6YTvkB66F7i+BDE7tYhCV/2bsZlmLSdGUYHcaKalSspmz+f6Kuvpvvvftfh10deeinQMiOm+dAhwkePcnJC4QukoHuJQRX5WJWJvTEp5z7Zy61N8o113o9rys+n+N55hJx/Pj3//FCXPtQu+cODNOzaTfr77xOUfLadIIU/koLuJfKjElnSZxyNgf6xtGzasRIsm7cQNtz752TXrFyFCgwk+bm/YwoOBk6fe3/cmebeHz8/MeJSnqvbytLrfs7d43/N/sevck1o4ZVkDN1LfJ18Af/K9J8x1N9ufZejTzxhdAyniPv5LfT5ZDnmXl3vUZeEx/PC0GsZWFnAtH1rnJBO+BIp6F7AWllJZJNvz8s+1bb486jfuRN7nfe+78b9+6nfsRPAqUs1fJl8Aet7DOTm3Z/QVFDgtOsK7ydDLl6g6t33WLz8Ga6f8mdqg8KMjuMW27v35frvV2HZvOXEzTbeRNvtlNz/AId27+OnV9zv3M28leIfw2bww7yvGBId7bzrCq8nBd0LWDZu5GBkgt8Uc4CcbmlgNmPZsN4rC/qxj5ZQv3Urr15wfYeK+ZnG1k9VHhrNwszp3B8T08mEwhfJkIuH01Yr9Zs2sT3ev9bzaAwMJjQzE8uGjUZH6TBbTQ1Hn36a0KFDWZk6wqVt1e/YwcE5c7A3NLi0HeEdpKB7uIZdu7BbLOyI72N0FLdL/MujpLz8ktExOqxs/gvYysvp8eCDaOXaHzF7bS21X6yk4tVXXdqO8A4y5OLhLBs2ALAjzr966AD9/7mr3ec9fVndgNhYYm+8kdDMwYBrP7QMv/BCIi6fQPlLLxNz/fUEduvm0vaEZ5OC7uEiJ00msEcPqr71zz+mrs5bjV0plvS92OgoZ3Xy2Hevln/ctOZOwty57F81jfJ/LqTHffPc0qbwTFLQPVxQcq+WOwK/9d0Fuc5mWOn3JNeWenxBB0itPkyv2jLWJZ4Pblri+PgvkjuTs7jszUVMLUmmNCzW4/+KEa7hn90+L9FUWEjVBx9iq601OophNif0p1ddGYm1ZUZHOadbcpYxd8s7hFnd/wHlogETeXXQFKqCI93etvAc0kP3MG3/dP/h919xa85SLvu6HkKiDExlnOwe/WEHjDiay9KIeKPjnNHA8nxGH9nNK4OmYDGHur390rBYPjxvvNvbFZ5FeugebFjp9xRGJFDpp8UcoDg8npKwOLKO7DE6yhlprfnpruVUBEeypM9YQ7NcdnATP9/5saEZhHGkoHsos62ZzLJ9bE7ob3QUYynF2sTzsZkCwEO3Yatbu5Yh5ftZ3H8CjYHBhmZJrTnCtXmrady/39AcwhhS0D3U+RX5BNutbE7IMDqK4V7KnMYjo3/qtg8aO0o3NbEjLp1Pe48xOgof9L2ExgAzZQsWGB1FGEAKuoc6r6qIZhXADj+7Q/RsgmzNRkdoV+Rll/H7i39Nc4DxH0lVB4ezNP0iqpcuo/GA723jJ85OCrqHei/jMm6c9EcaDP4T3lPckrOUBSuf9KhhF22zUfX++9gbG42OcpL/nHcJKiiI8n8uNDqKcDPjuxTijKqDw42O4DGKIhJItFTQ91ix0VFOqF66lJIHHsQU5VkfWleFRJJw112YkxKNjiLcTHroHujC4p08uP41v1sD/WzW9xyEDcVFJTuMjgK0jJuX/uN5QgYNInLCBKPjnKbbT24i8vLLjY4h3EwKugcac3gnQ8r2UWfAfGZPdSw4gl1x6VxYstPoKABUvf8+zUVFdJ97J8rkmT9GtpoaSv/xPE1FRUZHEW7i0P+JSqlJSqlcpVSeUuq0xSKUUgOUUuuUUo1KqbudH9N/KG1n5JE9bEroj93FK/V5m7WJg0mvPmz4Lj32+nrKXniR0KwRhI/z3LXa7RYL5QsXUr7wX0ZHEW5yzoqhlAoA5gOTgUHALKXUoFNOqwDuAJ5yekI/06/yILGNtWzoeeq3WKxJGsKLmVcTYPAuPdbycgITE0m4806Uh06lBDD36EHMdTOo+uADmos957MH4TqOdAFHAXla6/1a6yZgMTC97Qla66Na642AZ84r8yKjDu/GhmJjjwFGR/E4ZWExLOk7jgCDd+kJSk4m7Z3FhGVlGZrDEXG33gpA+UsvG5xEuIMjBb0XcLDN46LW5zpMKTVbKZWtlMouLS3tzCV8XkVoFJ/3HulX2811RLC1kar336cpP9+Q9mu++gprRYVH98zbMiclEXP11VS9+y7NR44YHUe4mCMFvb3/czs1GVhrvVBrnaW1zurevXtnLuHzlqVfxN8v+JHRMTxWiK2Jkj/+iar3/+P2tpuPHuXQnXM5+qR3jSzG3Tab8HHj0LJNnc9zZB56EZDS5nEyIANyLtB8+DBmm9Uj7jj0VMeCIwkfN5ZjS5e6fYZJ2Ysvoq1W4n95u9va7KzTNpuOvhJe3En+Y72NCSTcwpGfho1AhlIqXSkVBMwElrg2ln8q+eMfeebr54yO4fGir5qGtaTErRtINxUUUPXue8T+6DqCUlPd1q4z9awrp/rTz4yOIVzonF1BrbVVKTUH+AwIAF7RWucopW5vPb5AKdUTyAaiALtS6k5gkNa62nXRfYutupq6dd+xtfdFRkfxeJETfoApKoqq/3uH8DGj3dJm6d+fQ5nNTK7IoNJNW8s524/3rKB4zS7CRo2UvUd9lEN/r2qtl2ut+2mt+2qt/9L63AKt9YLWrw9rrZO11lFa65jWr6WYd0DNylXQ3MzqXkONjuLxTKGhxFxzDc0lh9E2m8vb01YrWtuJ+9lPvXpt+nf6TUA3NFDx6mtGRxEuIoO1HqL6008wJyWxNybl3CcLEn53FyooyC1tqcBAkp95Bm23w/2fuKVNVyiKTCBq8mQqFy2i2y0/IzA21uhIwsnkVkQPYDt2jLq164icNMlj1/z2NMeLua2qCm21uqyd2m/W0JiX19Kmh97i3xHxv7wdu8VCxRtvGB1FuID00D2AKSqKtLfeJCAuDp7fZnQcj3d8Bkf6sWKe+fo5nr3gR3yVMtzpO93bamoonjeP4PR0er/1plOvbZTgjAyipk5FNzYZHUW4gBR0D6CUInTo8bFzKeiOyo/qSUl4HNfvXcnXycOceu20ecv49bb3mVxewa8ybyLPSz8IbU/Sk094zY1RomO8/29IL2ctL6fkj3+iqbDQ6CheRysT7/SbQFrNES4syXHqtQeX7WPqgXV82Pdi8mKSnXptox0v5pYtW7DV1BicRjiT9NANcnzY4Oq81dy2cwnXl6VSGNXT4FTeZ3Wvody4ZwU37vkMbb0bFdj1/6XtDQ38dsu7lIR1482B/+OElJ6nqaCAglk3EHfbbSTMvdPoOMJJpIdusMsLs8mNSZFi3kl2UwCvDppC7+ojWDZtds5FlWJd0mCeGzaDRh/dAjCod2+irrySitdfp/noUaPjCCeRgm6gPlWH6FtdzOepnr9qnyf7NimT2ybcQ/joUV2+ltYaU3Awr5w/la0J/ZyQznN1/+0daKuVsvkvGB1FOIkUdANNLMym2RTA18kXGB3FuylFUWQCAI3ff9/pyzQVFnLg6muo3+EZuyK5WlBqKrE/+hFV771H44EDRscRTiBj6AaqM4fwRUqWLJXrJLXffMPBX8ym17PPEDVpUodea62s5ODs27BVVhIQ7b13g57LqYt2xTRk8PegSJIPHCA4Pd2gVMJZpKAb6C0f/cDNKOFjxhAyZAglDzxIUHofQvo7NmRiq6mh6PZf0nzoEKmvvdq6+JZnbEbtalUhkfxs4n3s/8EPjI4inECGXAygtWZART7oTi0rL84g/Q8ruC5xGkdtAWy8/ieMnfPa6cvInsJWXU3hT39G/a5dJP3tacJGjHBTWs9hNwWg7XaqV6xwy9o4wnWkoBugfssWnln9PJcWbTE6is8pC43hDxfeSqDdyhPfvECItfGs55vCwgiIjSX5H88RNXGim1J6ntrVqzl0x2+p+o/7Nw4RziMF3QAVr71OjTmUdYnnGx3FJ+VHJ3H3+Dks7j+BhtZph3UbNtB85CjNxcVUr1hB0W/uwFpaigoMJOVfC4m89FJjQxss4pJLCB0xgtK/PYPt2DGj44hOkjF0N2vcv5+azz9nacYPfHaOsyc4FNGdQxEt2xxaNm+m8Cc3n3Q8IC6Oht27Gfz0BiPieRylFD0ffIADP5xB6T+ep+eDDxgdSXSCFHQ3K1/4L1RwMB/1vdjoKH4jZOBAUl95uXV5BUXweX0JHTas5a7S5b6zRktXhQwcSOzM66l8+22ip11F6JAhRkcSHSQF3Y3sjY3UbVhPzI+u41h1hNFx/EafP69q/Sqy5V9bj8B7shVbe7rfdRcNuXvRTbIaozeSgu5GpuBg+n76KbqxEf6y2ug4QpwmICKCtEVvGR1DdJIUdDexlpVhiorCFBQEbtppRwhHnTq902xr5obcL1iTNISVC2YblEp0lMxycZOSBx6kYOYstMw9F14g2NbM5YUbmZf9Fva6OqPjCAdJQXeD2m+/pfbrr4maMlk2FhBeoTYojCeyfkxibRmHH37E6DjCQTLk4mL97/6AF1c9jS08nmm58TT70M43wrftiO/L2wMmctNHHxE2Zgwx11xtdCRxDtJDd7FZuV+QaCnnH8Nm0BxgNjqOEB2yuP/lhI0axdHHHsNWK0Mvnk566C6k7XYGVBawIjWL7d3PMzqOEB1mVyZ6PfsMzSUlBESEGx1HnIMUdBdSJhP3XzSbILvV6ChCdFpgt24EdusGwLGly4gYfzEBUb67xLA3k4LuAlprKl55lehrrsZuCqDBFGB0JCE67fiUxh515fzriyfYG5vCgxf9gj1PXWtwMnEqGUN3gbL5L3D0ySepXrrU6ChCOM2R8Dgez/oxAyoK+NN3r8h0Rg8kBd1J0uYtI23eMn5+w58pe/55Pk/JYnhON6NjCeFU3/YawtMjZpJZtp+Cn92CtbLS6EiiDSnoTnTZwc3M3fwOW+P78twFM0DmnAsf9GXKCB4dfTONubnUrV1rdBzRhoyhO0mg3crM3C/YEd+Hh8bcgtUk31rhu75LHEzfxz/F3LMnAM0lJZgTEw1OJaTqdJG2WtE2G1ZTIPeNvY06cyiNgbJWi/B9Gc9uAiDtWDF///o5VqaM4KXBU9n19AyDk/kvKehd0Lh/P8Xz7iOod28wjaciNNroSEK43aGIBJb0Gcs1eavJOrqH6k/CiZw0SZa5MICMoXeCrbaWo88+y4FrrqW5oIDIH1wm4+XCbzUHBPLy4Kv43fg51JrDODT3Lgpuukk2nDaA9NA7qG79Bg7ddRe28nKirryShHt/jzkhAb6WNVqEf8vt1ps5l81l82gr1rIyVEAAWmtqPltB+LhxcqepG0hBPwd7YyN1a9YQEBND2IgRBPVOJWTgQLrf8RvZokuIU9iViWEbgoAkmLeMfpWF/P3r52gyBRI34VIir/gfwi8cQ2B8vNFRfZJyZH1updQk4O9AAPCS1vqxU46r1uNTAAvwU6315rNdMysrS2dnZ3c2t0vVrV2LZdNmLNnZ1G/dim5sJGrKFHr97enTNgIQQpyF1gysKGD8oa3MqMnFWloKQOobrxM+ahRNBw9iPXqU4H79CIiMNDisd1BKbdJaZ7V77FwFXSkVAOwFJgJFwEZgltZ6V5tzpgC/oaWgjwb+rrUefbbrOqOga63BZmsZq7PZMIWFAS1j3HaLBex2tNWGbmxANzcTMmAAAHXr1tGUn4+1rBxraSnWo0cJiI0l6X//CsCBa39I3a7d7I9OYmd8H7J7DGBb/HnY5BZ+ITrNpO2cV1XE0NI8Pu4zlobAYG7Ys4Kb9qxoOR4VhTkxEXNSEklPPkFARASWTZtoys/HFBVFQGQkKjgYU3AwIYMGAWCrqQGbDQLNKJMCkwlMppadwWipEb724ezZCrojQy6jgDyt9f7Wiy0GpgO72pwzHXhDt/x2+E4pFaOUStRal3Qxe7sOP/IolW+/DW1+GangYAZs29py/OGHqV7y8UmvqQoKZ9aUPwPwh+9e5aLDOQBUBkdQERJFYWQPnmjtfSclXUVV+iws5lBXxBfCL9mVib2xqeyNTT3x3PK0C/k+JpneNUfobqkk4VgVcYf3Mv6RLznw+FUc+/hjqha/c9J1VFAQA7ZvA+DIo3/h2EcfnXQ8IDaWfutabngqmvMbaleuPFHolVKYU1Ppu6xlWY7C227Dsn7DSa8P7teP9P9rabPgxpuoz8k56XjosKH0fvVVoKXz13jgwEnHw8deRMrzzwOwb9Jkmo8ePel41MTLSXr8cQe/ax3jSA99BjBJa31r6+ObgNFa6zltzlkKPKa1XtP6eCVwr9Y6+5RrzQaOb1DYH8h11htxo3igzOgQbibv2ff52/sF733PvbXW3ds74EgPvb2/V079LeDIOWitFwILHWjTYymlss/0546vkvfs+/zt/YJvvmdH5qEXASltHicDxZ04RwghhAs5UtA3AhlKqXSlVBAwE1hyyjlLgJ+oFmOAY64aPxdCCNG+cw65aK2tSqk5wGe0TFt8RWudo5S6vfX4AmA5LTNc8miZtvgz10U2nFcPGXWSvGff52/vF3zwPTs0D10IIYTnk7VchBDCR0hBF0IIHyEFvQuUUncrpbRSyqcXplBKPamU2qOU2q6U+kApFWN0JldRSk1SSuUqpfKUUvOMzuNqSqkUpdSXSqndSqkcpdRvjc7kLkqpAKXUltb7aHyCFPROUkql0LIcQqHRWdzgc2Cw1noILctA3GdwHpdoXeZiPjAZGATMUkoNMjaVy1mB32mtBwJjgF/7wXs+7rfAbqNDOJMU9M57Bvg97dxA5Wu01iu01tbWh9/Rcp+BLzqxzIXWugk4vsyFz9JalxxfSE9rXUNLgetlbCrXU0olA1cCLxmdxZmkoHeCUmoacEhrvc3oLAa4BfjE6BAu0gs42OZxEX5Q3I5TSqUBFwDrDY7iDs/S0iGzG5zDqWQ99DNQSn0B9Gzn0APA/cAV7k3kWmd7v1rrj1rPeYCWP9EXuTObGzm0hIUvUkpFAO8Dd2qtq43O40pKqanAUa31JqXUpQbHcSop6Gegtb68veeVUplAOrCtdVnOZGCzUmqU1vqwGyM61Zne73FKqZuBqcAE7bs3L/jlEhZKKTMtxXyR1vo/Rudxg7HAtNZlv0OAKKXUW1rrGw3O1WVyY1EXKaXygSyttTeu2uaQ1g1O/gZcorUuNTqPqyilAmn50HcCcIiWZS9u0FrnnPWFXqx1c5rXgQqt9Z0Gx3G71h763VrrqQZHcQoZQxeOeB6IBD5XSm1VSi0wOpArtH7we3yZi93A//lyMW81FrgJ+EHrf9utrT1X4YWkhy6EED5CeuhCCOEjpKALIYSPkIIuhBA+Qgq6EEL4CCnoQgjhI6SgCyGEj5CCLoQQPuL/AQ/FGXux1HZjAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = 50\n",
    "x_min = -5.0\n",
    "x_max = 5.0\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.hist(samples, bins=bins, range=(x_min, x_max), density=True, color='tab:blue', label=\"samples\")\n",
    "xs = np.linspace(x_min, x_max, 1000)\n",
    "ys = inv_dist_fun(xs)\n",
    "ax.plot(xs, ys, color='tab:red', linestyle='--', label=\"target\")\n",
    "ax.set_ylim([0.0, 0.5])\n",
    "ax.legend(loc=\"upper right\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "02d3c230780613a338b307efe1c8520401457e802ba0c889aabda605c4f627d4"
  },
  "kernelspec": {
   "display_name": "Python 3.8.11 64-bit ('base': conda)",
   "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.8.11"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
	{
	"cells": [
	{
	"cell_type": "code",
	"execution_count": 1,
	"metadata": {},
	"outputs": [],
	"source": [
	"import numpy as np\n",
	"import matplotlib.pyplot as plt"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 2,
	"metadata": {},
	"outputs": [],
	"source": [
	"# Parameters\n",
	"n_mutate = 10000"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 3,
	"metadata": {},
	"outputs": [],
	"source": [
	"def inv_dist_fun(x):\n",
	" r = 0.3\n",
	" mu1 = 1.0\n",
	" mu2 = -2.0\n",
	" x1 = x - mu1\n",
	" x2 = x - mu2\n",
	" e1 = np.exp(-x1 * x1) / np.sqrt(np.pi)\n",
	" e2 = np.exp(-x2 * x2) / np.sqrt(np.pi)\n",
	" return r * e1 + (1.0 - r) * e2"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 4,
	"metadata": {},
	"outputs": [],
	"source": [
	"class HMC(object):\n",
	" def __init__(self, x, inv_dist_fun, sigma, seed=0):\n",
	" self.rng = np.random.RandomState(seed)\n",
	" self.x = x\n",
	" self.inv_dist_fun = inv_dist_fun\n",
	" self.sigma = sigma\n",
	" self.L = 50\n",
	" self.delta = 0.03\n",
	"\n",
	" def sample(self):\n",
	" x0 = self.x\n",
	" u0 = self.rng.normal(0.0, 1.0)\n",
	" eps = 1.0e-12\n",
	"\n",
	" # Leap-frog\n",
	" x1 = x0\n",
	" u1 = u0\n",
	" for l in range(self.L):\n",
	" x_temp = x1 + 0.5 * self.delta * u1\n",
	" x_plus = np.log(self.inv_dist_fun(x_temp + eps))\n",
	" x_minus = np.log(self.inv_dist_fun(x_temp - eps))\n",
	" dx = (x_plus - x_minus) / (2.0 * eps)\n",
	" u1 = u1 + self.delta * dx\n",
	" x1 = x_temp + 0.5 * self.delta * u1\n",
	"\n",
	" exp_H0_minus = self.inv_dist_fun(x0) * np.exp(-0.5 * u0 * u0)\n",
	" exp_H1_minus = self.inv_dist_fun(x1) * np.exp(-0.5 * u1 * u1)\n",
	" a = min(1.0, exp_H1_minus / (exp_H0_minus + 1.0e-12))\n",
	" if self.rng.random() < a:\n",
	" self.x = x1\n",
	" else:\n",
	" self.x = x0\n",
	"\n",
	" return self.x"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 5,
	"metadata": {},
	"outputs": [],
	"source": [
	"sampler = HMC(0.0, inv_dist_fun, 1.0)\n",
	"\n",
	"# Burn-in\n",
	"n_burnin = max(n_mutate // 100, 100)\n",
	"for i in range(n_burnin):\n",
	" sampler.sample()\n",
	"\n",
	"# Sampling\n",
	"samples = []\n",
	"for i in range(n_mutate):\n",
	" samples.append(sampler.sample())"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 6,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAulUlEQVR4nO3deXxU1d3H8c+ZZLKvJAQSkpCAYZMAQlgURC3iA4igFitYra21aFtqxWrFpa1V28e1WitKqbti8VGrIqCioCKCQNgJEAyQhJAAWcky2WbmPH8k0AABJsnM3Fl+79fLl5m5d+75TjS/nJw59xyltUYIIYT3MxkdQAghhHNIQRdCCB8hBV0IIXyEFHQhhPARUtCFEMJHSEEXQggf4VBBV0pNUkrlKqXylFLz2jl+qVLqmFJqa+s/f3R+VCGEEGcTeK4TlFIBwHxgIlAEbFRKLdFa7zrl1G+01lNdkFEIIYQDHOmhjwLytNb7tdZNwGJgumtjCSGE6Khz9tCBXsDBNo+LgNHtnHehUmobUAzcrbXOOfUEpdRsYDZAeHj4iAEDBnQ8sRBC+LFNmzaVaa27t3fMkYKu2nnu1PUCNgO9tda1SqkpwIdAxmkv0nohsBAgKytLZ2dnO9C8EEKI45RSBWc65siQSxGQ0uZxMi298BO01tVa69rWr5cDZqVUfCeyCiGE6CRHCvpGIEMpla6UCgJmAkvanqCU6qmUUq1fj2q9brmzwwohhDizcw65aK2tSqk5wGdAAPCK1jpHKXV76/EFwAzgl0opK1APzNSyjKMQQriVMqruyhi6EP6lubmZoqIiGhoajI7iFUJCQkhOTsZsNp/0vFJqk9Y6q73XOPKhqBBCdFlRURGRkZGkpaXROkIrzkBrTXl5OUVFRaSnpzv8Orn1XwjhFg0NDcTFxUkxd4BSiri4uA7/NSMFXQjhNlLMHdeZ75UUdCGE8BFS0IUQfqGqqooXXnjB5e18+OGH7Np16lJX7iEFXQjhFzpa0LXW2O32DrcjBV0IIVxs3rx57Nu3j2HDhjF37lwmTJjA8OHDyczM5KOPPgIgPz+fgQMH8qtf/Yrhw4dz8OBBHnnkEQYMGMDEiROZNWsWTz31FAD79u1j0qRJjBgxgosvvpg9e/awdu1alixZwj333MOwYcPYt2+fW9+jTFsUQhii4KafnPZc5ORJdLvhBuz19Rycfdtpx6OvuYaYa6/BWlnJoTt+e9Kx3m++cdb2HnvsMXbu3MnWrVuxWq1YLBaioqIoKytjzJgxTJs2DYDc3FxeffVVXnjhBbKzs3n//ffZsmULVquV4cOHM2LECABmz57NggULyMjIYP369fzqV79i1apVTJs2jalTpzJjxozOfms6TQq6EMLvaK25//77Wb16NSaTiUOHDnHkyBEAevfuzZgxYwBYs2YN06dPJzQ0FICrrroKgNraWtauXct111134pqNjY1ufhenk4IuhDDE2XrUptDQsx4PjI09Z4/8bBYtWkRpaSmbNm3CbDaTlpZ2Ys53eHj4ifPOdCe93W4nJiaGrVu3djqDK8gYuhDCL0RGRlJTUwPAsWPHSEhIwGw28+WXX1JQ0P6KtOPGjePjjz+moaGB2tpali1bBkBUVBTp6em8++67QEvh37Zt22ntuJsUdCGEX4iLi2Ps2LEMHjyYrVu3kp2dTVZWFosWLeJMm+2MHDmSadOmMXToUK699lqysrKIjo4GWnr5L7/8MkOHDuX8888/8cHqzJkzefLJJ7ngggvc/qGoLM4lhHCL3bt3M3DgQKNjdFhtbS0RERFYLBbGjx/PwoULGT58uFvabu97JotzCSFEJ82ePZtdu3bR0NDAzTff7LZi3hlS0IUQ4izefvttoyM4TMbQhRDCR0hBF0IIHyEFXQghfIQUdCGE8BHyoagQwhBp85Y59Xr5j13p1Ot1xaWXXspTTz1FVla7swtdRnroQgjhI6SgCyH8Ql1dHVdeeSVDhw5l8ODBvPPOOzz88MOMHDmSwYMHM3v27BNrt1x66aXMnTuX8ePHM3DgQDZu3Mi1115LRkYGDz74INCy1O6AAQO4+eabGTJkCDNmzMBisZzW7ooVK7jwwgsZPnw41113HbW1tUDLcr6DBg1iyJAh3H333U55j1LQhRB+4dNPPyUpKYlt27axc+dOJk2axJw5c9i4cSM7d+6kvr6epUuXnjg/KCiI1atXc/vttzN9+nTmz5/Pzp07ee211ygvLwdaltqdPXs227dvJyoq6rQNNMrKynj00Uf54osv2Lx5M1lZWfztb3+joqKCDz74gJycHLZv337il0RXSUEXQviFzMxMvvjiC+69916++eYboqOj+fLLLxk9ejSZmZmsWrWKnJycE+cfXx89MzOT888/n8TERIKDg+nTpw8HDx4EICUlhbFjxwJw4403smbNmpPa/O6779i1axdjx45l2LBhvP766xQUFBAVFUVISAi33nor//nPfwgLC3PKe5QPRYUQfqFfv35s2rSJ5cuXc99993HFFVcwf/58srOzSUlJ4aGHHjqxhC5AcHAwACaT6cTXxx9brVYAlFIntXHqY601EydO5N///vdpeTZs2MDKlStZvHgxzz//PKtWrerye5QeuhDCLxQXFxMWFsaNN97I3XffzebNmwGIj4+ntraW9957r8PXLCwsZN26dQD8+9//Zty4cScdHzNmDN9++y15eXkAWCwW9u7dS21tLceOHWPKlCk8++yzTltXXXroQghDuHua4Y4dO7jnnnswmUyYzWZefPFFPvzwQzIzM0lLS2PkyJEdvubAgQN5/fXXue2228jIyOCXv/zlSce7d+/Oa6+9xqxZs07saPToo48SGRnJ9OnTaWhoQGvNM88845T3KMvnCiHcwluXzz2T/Px8pk6dys6dO13WRkeXz5UhFyGE8BFS0IUQohPS0tJc2jvvDCnoQgi3MWqI1xt15nslBV0I4RYhISGUl5dLUXeA1pry8nJCQkI69DqZ5SKEcIvk5GSKioooLS01OopXCAkJITk5uUOvkYIuhHALs9lMenq60TF8mgy5CCGEj3CooCulJimlcpVSeUqpeWc5b6RSyqaUmuG8iEIIIRxxzoKulAoA5gOTgUHALKXUoDOc9zjwmbNDCiGEODdHxtBHAXla6/0ASqnFwHRg1ynn/QZ4H+j4/bPCI5xpBxlP2glGCHFmjgy59AIOtnlc1PrcCUqpXsA1wIKzXUgpNVspla2UypZPuoUQwrkc6aGrdp47dSLps8C9WmvbqctHnvQirRcCC6FlLRcHMwon68hejqHNDZT8+c/ETJ9O6LBhrgslhOgyR3roRUBKm8fJQPEp52QBi5VS+cAM4AWl1NXOCCjcr0/VIS4q3gFAYl051R8tIX/mLI4+/TTabjc4nRDiTBzpoW8EMpRS6cAhYCZwQ9sTtNYnJpcqpV4DlmqtP3ReTOEuibVl/HXtQmqCwljfcxD7Y3qR8c1qjjz+BOX/eglMASTMvdPomEKIdpyzoGutrUqpObTMXgkAXtFa5yilbm89ftZxc+E9Auw27tv4JqD545ifYzMFAGAKD6fnnx8CrSn/5z8JHTaUyMsuMzSrEOJ0Dt0pqrVeDiw/5bl2C7nW+qddjyWM8MO8r8g4doi/jLyJkoj4k44ppej5h5aNbIP79DEinhDiHOTWfwFAdGMNM3NXsjZxMGt6DW33HBUUROIjD7s5mRDCUXLrvwCgJiic+UOv5ZXzzz3nvKnoEMX3zsNaUeGGZEIIR0kPXQBgVyZWpra7q9Vp0xxTao6wYOUSzL2S6H7HHe6IJ4RwgBR0wf/kryeqqY53My6Ds9xHcNzByB5s7DGA/i+9wUUHU2kOMJ90XO4sFcIYMuTi50zazqzczxl+NNehYn7ch+eNJ6apjkuLtrgwnRCiI6Sg+7lRh3fRo76KpeljO/S6rfHnURiZwP8UbHBRMiFER8mQi5+7omAj5SFRrEs8v2MvVIoP+1xMSu1RTNqOXUnfQAijSUH3YxFNFrKO7GFJn7HYW28i6ohP0i90QSohRGdJt8qPRTZb2JzQj69Shnf6GkrbGVy2D6VljRchjCYF3Y+VhMfz0IU/Jy+mYxvRtjX+0DaeXPMi/SoPnvtkIYRLSUH3U2HNDcTXV3X5OpsS+mNVJi4q2dn1UEKILpGC7qfGFW/nzc8eJbnmaJeuUxsUxvb481oKupYl7oUwkhR0PzW6JIejoTEURXTv8rXWJg0mubaU1JojTkgmhOgsKeh+yGxrZnjpXtb3HNShm4nOZF3PlimPF5bkdPlaQojOk2mLfmhYaR4htma+69nBuednUBEazW8vuYP90UlOuZ4QonOkoPuh0Yd3YQkMZkd8X6ddc29sqtOuJYToHCnofmjRgIl8m5RJc4Dz/vOHNjcwa+8XbOneD5DFuYQwghR0P1QZEkVlSJRTr9kYYGZS/nqiGi1Ova4QwnHyoaifqf3mG6bu/xaTk+/stJsC2NI9gxFH96Bl+qIQhpCC7mcq33mHH+Z9hZ2uz2451aaE/sQ3VNP4/fdOv7YQ4tykoPsRbbNhWb+Brd0znDJd8VSbe/QHoG7Nt06/thDi3KSg+5GGnBzsNTUtBd0FykJj2B3bG93U5JLrCyHOTj4U9SN1674DYGv381zWxl2X/Ib822WWixBGkB66H2kuKSZ44ECOBUe6vC35YFQI95OC7kcSH3qI9HcWu7SNYGsj+6ZcScVrr7u0HSHE6aSg+xkVFOTS6zcGBoPdjmWD7DUqhLtJQfcTFW+8SeGtv0A3N7u8rbBRo7BkZ6NtNpe3JYT4LynofqJ2zTc0Hy5Bmc0ubyts1CjsNTU07Nnj8raEEP8lBd0PaJuN+i1bCRs+wi3thY0cCYBlvQy7COFOUtD9QGNeHvaaGsKy3FPQzT0S6PbzWwju388t7QkhWsg8dB+VNm/Zia+n7v+WXwMTPz/G0bXLzvwiJ+pxzz1uaUcI8V/SQ/cDVcERfJM0hKNhsW5rU2tN44ED2I4dc1ubQvg76aH7gTW9hrKm11C3tZc2bxm9q0tYsOppnho+k5WpWQDkPyZ3kArhStJD93HB1iZCrI1ub7cwsgd1gSEMrMh3e9tC+CuHCrpSapJSKlcplaeUmtfO8elKqe1Kqa1KqWyl1DjnRxWdMbZ4B+8t+wNJtaVubVcrE7u79WZQRYFb2xXCn52zoCulAoD5wGRgEDBLKTXolNNWAkO11sOAW4CXnJxTdNLAynwaA8wcDo9ze9u7u/Wmd/Vhwprr3d62EP7IkR76KCBPa71fa90ELAamtz1Ba12r/7saUzggKzN5iH6VB9kbk4JduX90bVdcOiY0AyoL3d62EP7IkZ/yXsDBNo+LWp87iVLqGqXUHmAZLb300yilZrcOyWSXlrp3CMAfmW3N9DlWTG5sqiHt745N5ZFRN7M3JsWQ9oXwN44U9Pa2tjmtB661/kBrPQC4GnikvQtprRdqrbO01lndu3fvUFDRcX2PFROo7eyNNaagNgYGszYpk9qgMEPaF8LfOFLQi4C2FSEZKD7TyVrr1UBfpVR8F7OJLioNjeGfg6eRE5dmWIbE2jKuyfva6ZtSCyFO50hB3whkKKXSlVJBwExgSdsTlFLnKdWySaVSajgQBJQ7O6zomPLQaD48b7xbNrQ4kwGVhcze+TG9qw8blkEIf3HOG4u01lal1BzgMyAAeEVrnaOUur31+ALgh8BPlFLNQD1wvZYtaww38vBu8mJ6URkSZViGXd16AzBI5qML4XIO3SmqtV4OLD/luQVtvn4ceNy50URXRDRZePi7l3l10GT+r98Ew3IcCetGVVA4/SoPnvtkIUSXyJ2iPqpfVUsB3WvQDJcTlCI3NpX+MnVRCJeTgu6jjveI98YkG5wE9samkFRXhr1ebjASwpVkcS4f1b+ykMKIBCzmUKOj8GHfi3kv4zL2hhqfRQhfJj10H6S1brlD1KD556eymENpCnD91ndC+Dvpofuou8f/GuVBE42uzltN6XP76X7Hb4yOIoTPkh66D1JKURIeT3GE59yN26/qIFXvvWd0DCF8mhR0H1T9ySdckb/e6Bgn2RObivXoUZqPHDE6ihA+Swq6D6p8+99MLvjO6BgnOT59sn77doOTCOG7pKD7GG2zUZ+TQ25sb6OjnGRfdBKYzTRIQRfCZaSg+5jGvH1oi4VcD5nhclxzgJnwUaOMjiGET5NZLj6mYUdLD9jwO0TbkfqybGQlhCtJD93HNBUVYYqOptiALeeEEMaSgu5jEu68k4yvvkQbsOXcuVjLyth/1VVUffCh0VGE8Eme91MvuszkobfYB3TrRvPhI9Rv3Wp0FCF8khR0H1K/YyeFt91G44EDRkdplzKZCM0cTP0OmekihCtIQfch9Zs3Uff1akzh4UZHOaOQzCE05u6VlReFcAGZ5eJD6rfvIDAxEXNCgtFR2pU2bxmjS5p4yGbjqjteJicuHYD8x640OJkQvkEKupdLm7fsxNevfPkd+2KSmNjmOU+TG5vKyuTh1AcEGR1FCJ8jQy4+IqqxjkRLObkeOP+8raqQSJ7KuoH9Mb2MjiKEz5GC7iMim+rYEdeH3d3SjI5yblrT3VJpdAohfI4UdB9xKDKB31/8qxPj0p5s2v41vLHiL0Q31hodRQifIgXdR5i03egIDtsX3TLc0k82jhbCqaSg+wKtef2zR7lhzwqjkzgkL6YXNhQDpKAL4VRS0H1AD0sF8Q3VVAVHGB3FIY2BwRRE9aRf5UGjowjhU6Sg+4D+rYXR02e4tJUbm0r/ykLwoH1PhfB2Mg/dB/SrLKTRFEh+VKLRURz2We9RbEnI8KqxfyE8nRR0H9C/6iB5McnYTAFGR3FYbrfe5OJZuyoJ4e2koPuAb5KGYgkMNjpGh51XVUSotdHoGEL4DCnoPmBJ33FGR+iUX+xYQpC9GbjD6ChC+AT5UNTLdbdUEtNQY3SMTsmNTaXvsWLsTU1GRxHCJ0hB93I35H7OP1c+6ZWzRXJjUzHbbTTu2WN0FCF8ghR0L9e/8iB7Y1NAKaOjdNjxaZb122TDCyGcQQq6F7PX1ZFafdir5p+3VRYaTXlIFPXbpaAL4QxS0L1Y/c4cAtBeW9BRigcu+gU9//Qno5MI4ROkoHux+u3bAO+6Q/RUBVGJBER47pZ5QngThwq6UmqSUipXKZWnlJrXzvEfK6W2t/6zVik11PlRxamip0zhryNvpDrYewtiZFMdpc89R/2OnUZHEcLrnXMeulIqAJgPTASKgI1KqSVa611tTjsAXKK1rlRKTQYWAqNdEVj8l7lXL77pNczoGF1iUybKXlwAgYGEZg42Oo4QXs2RHvooIE9rvV9r3QQsBqa3PUFrvVZrfXwLmu+AZOfGFKeyVlZS9d57RDd65xz04yzmUIL69KFh+w6jowjh9Rwp6L2AtuucFrU+dyY/Bz5p74BSarZSKlsplV1aWup4SnEaS3Y2JQ/+gcS6CqOjdFnokCHUb9+O9sK59EJ4EkcKensTnNv9yVNKXUZLQb+3veNa64Va6yytdVb37t0dTylO07B9O5jN7ItOMjpKl4UOycRWUUHzoWKjowjh1Rwp6EVASpvHycBpP3lKqSHAS8B0rXW5c+KJM6nftp2QAQNoDjAbHaXLQoYMwRQWRnORbHghRFc4UtA3AhlKqXSlVBAwE1jS9gSlVCrwH+AmrfVe58cUbWmbjYadOwnNzDQ6ilOEDBxIv40bCB8zxugoQni1c85y0VpblVJzgM+AAOAVrXWOUur21uMLgD8CccALquUWdKvWOst1sf1bU0EBdouF0KFDYJ3RabpOmeR2CCGcwaHlc7XWy4Hlpzy3oM3XtwK3OjeaOJPgPn3IWPstKigY1n1ldBynqP78c8oX/ovei97CFBRkdBwhvJJ0jbxUYLduPneHZcOOHTTszDE6hhBeSwq6Fzr8179S/Um7M0O9VtgFFwBQv2WzwUmE8F5S0L2M3WKh8q1FNH7/vdFRnCowPh5z71Qsm7cYHUUIryUF3cs05OSA3U7IkCFGR3G6sAuGU795s9xgJEQnyZ6iXqZ+W8sKi6FDfW/9s4hLxqObGrHXWXzu8wEh3EEKupexbN5CUFoagbGxRkdxuqjJk4maPNnoGEJ4LSno3kZrwsb41kKWafOWnfQ4rLkeizmU/MeuNCiRezUVHaJu7bc0FxYSMWHCiQ+IhegoKeheJuXFF3x6jPl3m/5N/8pCZl/e7nJAPqWpoIAjTzxJ7apVoDXKbCZsdMvdspbNmyl/+RUS7ppLcN++BicV3kI+FPVCygs3hHbUwcgEUmpLiWqsMzqKS9lq6zhw3Y+wrF9P/C9vp8/y5fTfvo2Ii8cB0HyoGMuGDRy45loq3lrk07/EhfNID92LHH74EZqLi0lZ8KLRUVxmV7c0AAZW5Buaw9UCIsJJfPQRQocOw9wj4bTj0VdNJfzCMZQ8+AeOPPoozQcLSbj3XlkmQZyVFHQvkTZvGf9c+QUlYXE8dMqYsy/ZG5tKswrg/IoDRkdxOq01R/73fwkdMpToqVcSdcUVZz0/MD6e5Bfmc/Txx6l4/Q2CBw4k5uqr3RNWeCX5de8lIpospNYcZVdcmtFRXKopwMy+mF4MKs83OorTlT0/n8o33qRxz26HX6NMJnrcdx/JC14keto0F6YTvkB66F7i+BDE7tYhCV/2bsZlmLSdGUYHcaKalSspmz+f6Kuvpvvvftfh10deeinQMiOm+dAhwkePcnJC4QukoHuJQRX5WJWJvTEp5z7Zy61N8o113o9rys+n+N55hJx/Pj3//FCXPtQu+cODNOzaTfr77xOUfLadIIU/koLuJfKjElnSZxyNgf6xtGzasRIsm7cQNtz752TXrFyFCgwk+bm/YwoOBk6fe3/cmebeHz8/MeJSnqvbytLrfs7d43/N/sevck1o4ZVkDN1LfJ18Af/K9J8x1N9ufZejTzxhdAyniPv5LfT5ZDnmXl3vUZeEx/PC0GsZWFnAtH1rnJBO+BIp6F7AWllJZJNvz8s+1bb486jfuRN7nfe+78b9+6nfsRPAqUs1fJl8Aet7DOTm3Z/QVFDgtOsK7ydDLl6g6t33WLz8Ga6f8mdqg8KMjuMW27v35frvV2HZvOXEzTbeRNvtlNz/AId27+OnV9zv3M28leIfw2bww7yvGBId7bzrCq8nBd0LWDZu5GBkgt8Uc4CcbmlgNmPZsN4rC/qxj5ZQv3Urr15wfYeK+ZnG1k9VHhrNwszp3B8T08mEwhfJkIuH01Yr9Zs2sT3ev9bzaAwMJjQzE8uGjUZH6TBbTQ1Hn36a0KFDWZk6wqVt1e/YwcE5c7A3NLi0HeEdpKB7uIZdu7BbLOyI72N0FLdL/MujpLz8ktExOqxs/gvYysvp8eCDaOXaHzF7bS21X6yk4tVXXdqO8A4y5OLhLBs2ALAjzr966AD9/7mr3ec9fVndgNhYYm+8kdDMwYBrP7QMv/BCIi6fQPlLLxNz/fUEduvm0vaEZ5OC7uEiJ00msEcPqr71zz+mrs5bjV0plvS92OgoZ3Xy2Hevln/ctOZOwty57F81jfJ/LqTHffPc0qbwTFLQPVxQcq+WOwK/9d0Fuc5mWOn3JNeWenxBB0itPkyv2jLWJZ4Pblri+PgvkjuTs7jszUVMLUmmNCzW4/+KEa7hn90+L9FUWEjVBx9iq601OophNif0p1ddGYm1ZUZHOadbcpYxd8s7hFnd/wHlogETeXXQFKqCI93etvAc0kP3MG3/dP/h919xa85SLvu6HkKiDExlnOwe/WEHjDiay9KIeKPjnNHA8nxGH9nNK4OmYDGHur390rBYPjxvvNvbFZ5FeugebFjp9xRGJFDpp8UcoDg8npKwOLKO7DE6yhlprfnpruVUBEeypM9YQ7NcdnATP9/5saEZhHGkoHsos62ZzLJ9bE7ob3QUYynF2sTzsZkCwEO3Yatbu5Yh5ftZ3H8CjYHBhmZJrTnCtXmrady/39AcwhhS0D3U+RX5BNutbE7IMDqK4V7KnMYjo3/qtg8aO0o3NbEjLp1Pe48xOgof9L2ExgAzZQsWGB1FGEAKuoc6r6qIZhXADj+7Q/RsgmzNRkdoV+Rll/H7i39Nc4DxH0lVB4ezNP0iqpcuo/GA723jJ85OCrqHei/jMm6c9EcaDP4T3lPckrOUBSuf9KhhF22zUfX++9gbG42OcpL/nHcJKiiI8n8uNDqKcDPjuxTijKqDw42O4DGKIhJItFTQ91ix0VFOqF66lJIHHsQU5VkfWleFRJJw112YkxKNjiLcTHroHujC4p08uP41v1sD/WzW9xyEDcVFJTuMjgK0jJuX/uN5QgYNInLCBKPjnKbbT24i8vLLjY4h3EwKugcac3gnQ8r2UWfAfGZPdSw4gl1x6VxYstPoKABUvf8+zUVFdJ97J8rkmT9GtpoaSv/xPE1FRUZHEW7i0P+JSqlJSqlcpVSeUuq0xSKUUgOUUuuUUo1KqbudH9N/KG1n5JE9bEroj93FK/V5m7WJg0mvPmz4Lj32+nrKXniR0KwRhI/z3LXa7RYL5QsXUr7wX0ZHEW5yzoqhlAoA5gOTgUHALKXUoFNOqwDuAJ5yekI/06/yILGNtWzoeeq3WKxJGsKLmVcTYPAuPdbycgITE0m4806Uh06lBDD36EHMdTOo+uADmos957MH4TqOdAFHAXla6/1a6yZgMTC97Qla66Na642AZ84r8yKjDu/GhmJjjwFGR/E4ZWExLOk7jgCDd+kJSk4m7Z3FhGVlGZrDEXG33gpA+UsvG5xEuIMjBb0XcLDN46LW5zpMKTVbKZWtlMouLS3tzCV8XkVoFJ/3HulX2811RLC1kar336cpP9+Q9mu++gprRYVH98zbMiclEXP11VS9+y7NR44YHUe4mCMFvb3/czs1GVhrvVBrnaW1zurevXtnLuHzlqVfxN8v+JHRMTxWiK2Jkj/+iar3/+P2tpuPHuXQnXM5+qR3jSzG3Tab8HHj0LJNnc9zZB56EZDS5nEyIANyLtB8+DBmm9Uj7jj0VMeCIwkfN5ZjS5e6fYZJ2Ysvoq1W4n95u9va7KzTNpuOvhJe3En+Y72NCSTcwpGfho1AhlIqXSkVBMwElrg2ln8q+eMfeebr54yO4fGir5qGtaTErRtINxUUUPXue8T+6DqCUlPd1q4z9awrp/rTz4yOIVzonF1BrbVVKTUH+AwIAF7RWucopW5vPb5AKdUTyAaiALtS6k5gkNa62nXRfYutupq6dd+xtfdFRkfxeJETfoApKoqq/3uH8DGj3dJm6d+fQ5nNTK7IoNJNW8s524/3rKB4zS7CRo2UvUd9lEN/r2qtl2ut+2mt+2qt/9L63AKt9YLWrw9rrZO11lFa65jWr6WYd0DNylXQ3MzqXkONjuLxTKGhxFxzDc0lh9E2m8vb01YrWtuJ+9lPvXpt+nf6TUA3NFDx6mtGRxEuIoO1HqL6008wJyWxNybl3CcLEn53FyooyC1tqcBAkp95Bm23w/2fuKVNVyiKTCBq8mQqFy2i2y0/IzA21uhIwsnkVkQPYDt2jLq164icNMlj1/z2NMeLua2qCm21uqyd2m/W0JiX19Kmh97i3xHxv7wdu8VCxRtvGB1FuID00D2AKSqKtLfeJCAuDp7fZnQcj3d8Bkf6sWKe+fo5nr3gR3yVMtzpO93bamoonjeP4PR0er/1plOvbZTgjAyipk5FNzYZHUW4gBR0D6CUInTo8bFzKeiOyo/qSUl4HNfvXcnXycOceu20ecv49bb3mVxewa8ybyLPSz8IbU/Sk094zY1RomO8/29IL2ctL6fkj3+iqbDQ6CheRysT7/SbQFrNES4syXHqtQeX7WPqgXV82Pdi8mKSnXptox0v5pYtW7DV1BicRjiT9NANcnzY4Oq81dy2cwnXl6VSGNXT4FTeZ3Wvody4ZwU37vkMbb0bFdj1/6XtDQ38dsu7lIR1482B/+OElJ6nqaCAglk3EHfbbSTMvdPoOMJJpIdusMsLs8mNSZFi3kl2UwCvDppC7+ojWDZtds5FlWJd0mCeGzaDRh/dAjCod2+irrySitdfp/noUaPjCCeRgm6gPlWH6FtdzOepnr9qnyf7NimT2ybcQ/joUV2+ltYaU3Awr5w/la0J/ZyQznN1/+0daKuVsvkvGB1FOIkUdANNLMym2RTA18kXGB3FuylFUWQCAI3ff9/pyzQVFnLg6muo3+EZuyK5WlBqKrE/+hFV771H44EDRscRTiBj6AaqM4fwRUqWLJXrJLXffMPBX8ym17PPEDVpUodea62s5ODs27BVVhIQ7b13g57LqYt2xTRk8PegSJIPHCA4Pd2gVMJZpKAb6C0f/cDNKOFjxhAyZAglDzxIUHofQvo7NmRiq6mh6PZf0nzoEKmvvdq6+JZnbEbtalUhkfxs4n3s/8EPjI4inECGXAygtWZART7oTi0rL84g/Q8ruC5xGkdtAWy8/ieMnfPa6cvInsJWXU3hT39G/a5dJP3tacJGjHBTWs9hNwWg7XaqV6xwy9o4wnWkoBugfssWnln9PJcWbTE6is8pC43hDxfeSqDdyhPfvECItfGs55vCwgiIjSX5H88RNXGim1J6ntrVqzl0x2+p+o/7Nw4RziMF3QAVr71OjTmUdYnnGx3FJ+VHJ3H3+Dks7j+BhtZph3UbNtB85CjNxcVUr1hB0W/uwFpaigoMJOVfC4m89FJjQxss4pJLCB0xgtK/PYPt2DGj44hOkjF0N2vcv5+azz9nacYPfHaOsyc4FNGdQxEt2xxaNm+m8Cc3n3Q8IC6Oht27Gfz0BiPieRylFD0ffIADP5xB6T+ep+eDDxgdSXSCFHQ3K1/4L1RwMB/1vdjoKH4jZOBAUl95uXV5BUXweX0JHTas5a7S5b6zRktXhQwcSOzM66l8+22ip11F6JAhRkcSHSQF3Y3sjY3UbVhPzI+u41h1hNFx/EafP69q/Sqy5V9bj8B7shVbe7rfdRcNuXvRTbIaozeSgu5GpuBg+n76KbqxEf6y2ug4QpwmICKCtEVvGR1DdJIUdDexlpVhiorCFBQEbtppRwhHnTq902xr5obcL1iTNISVC2YblEp0lMxycZOSBx6kYOYstMw9F14g2NbM5YUbmZf9Fva6OqPjCAdJQXeD2m+/pfbrr4maMlk2FhBeoTYojCeyfkxibRmHH37E6DjCQTLk4mL97/6AF1c9jS08nmm58TT70M43wrftiO/L2wMmctNHHxE2Zgwx11xtdCRxDtJDd7FZuV+QaCnnH8Nm0BxgNjqOEB2yuP/lhI0axdHHHsNWK0Mvnk566C6k7XYGVBawIjWL7d3PMzqOEB1mVyZ6PfsMzSUlBESEGx1HnIMUdBdSJhP3XzSbILvV6ChCdFpgt24EdusGwLGly4gYfzEBUb67xLA3k4LuAlprKl55lehrrsZuCqDBFGB0JCE67fiUxh515fzriyfYG5vCgxf9gj1PXWtwMnEqGUN3gbL5L3D0ySepXrrU6ChCOM2R8Dgez/oxAyoK+NN3r8h0Rg8kBd1J0uYtI23eMn5+w58pe/55Pk/JYnhON6NjCeFU3/YawtMjZpJZtp+Cn92CtbLS6EiiDSnoTnTZwc3M3fwOW+P78twFM0DmnAsf9GXKCB4dfTONubnUrV1rdBzRhoyhO0mg3crM3C/YEd+Hh8bcgtUk31rhu75LHEzfxz/F3LMnAM0lJZgTEw1OJaTqdJG2WtE2G1ZTIPeNvY06cyiNgbJWi/B9Gc9uAiDtWDF///o5VqaM4KXBU9n19AyDk/kvKehd0Lh/P8Xz7iOod28wjaciNNroSEK43aGIBJb0Gcs1eavJOrqH6k/CiZw0SZa5MICMoXeCrbaWo88+y4FrrqW5oIDIH1wm4+XCbzUHBPLy4Kv43fg51JrDODT3Lgpuukk2nDaA9NA7qG79Bg7ddRe28nKirryShHt/jzkhAb6WNVqEf8vt1ps5l81l82gr1rIyVEAAWmtqPltB+LhxcqepG0hBPwd7YyN1a9YQEBND2IgRBPVOJWTgQLrf8RvZokuIU9iViWEbgoAkmLeMfpWF/P3r52gyBRI34VIir/gfwi8cQ2B8vNFRfZJyZH1updQk4O9AAPCS1vqxU46r1uNTAAvwU6315rNdMysrS2dnZ3c2t0vVrV2LZdNmLNnZ1G/dim5sJGrKFHr97enTNgIQQpyF1gysKGD8oa3MqMnFWloKQOobrxM+ahRNBw9iPXqU4H79CIiMNDisd1BKbdJaZ7V77FwFXSkVAOwFJgJFwEZgltZ6V5tzpgC/oaWgjwb+rrUefbbrOqOga63BZmsZq7PZMIWFAS1j3HaLBex2tNWGbmxANzcTMmAAAHXr1tGUn4+1rBxraSnWo0cJiI0l6X//CsCBa39I3a7d7I9OYmd8H7J7DGBb/HnY5BZ+ITrNpO2cV1XE0NI8Pu4zlobAYG7Ys4Kb9qxoOR4VhTkxEXNSEklPPkFARASWTZtoys/HFBVFQGQkKjgYU3AwIYMGAWCrqQGbDQLNKJMCkwlMppadwWipEb724ezZCrojQy6jgDyt9f7Wiy0GpgO72pwzHXhDt/x2+E4pFaOUStRal3Qxe7sOP/IolW+/DW1+GangYAZs29py/OGHqV7y8UmvqQoKZ9aUPwPwh+9e5aLDOQBUBkdQERJFYWQPnmjtfSclXUVV+iws5lBXxBfCL9mVib2xqeyNTT3x3PK0C/k+JpneNUfobqkk4VgVcYf3Mv6RLznw+FUc+/hjqha/c9J1VFAQA7ZvA+DIo3/h2EcfnXQ8IDaWfutabngqmvMbaleuPFHolVKYU1Ppu6xlWY7C227Dsn7DSa8P7teP9P9rabPgxpuoz8k56XjosKH0fvVVoKXz13jgwEnHw8deRMrzzwOwb9Jkmo8ePel41MTLSXr8cQe/ax3jSA99BjBJa31r6+ObgNFa6zltzlkKPKa1XtP6eCVwr9Y6+5RrzQaOb1DYH8h11htxo3igzOgQbibv2ff52/sF733PvbXW3ds74EgPvb2/V079LeDIOWitFwILHWjTYymlss/0546vkvfs+/zt/YJvvmdH5qEXASltHicDxZ04RwghhAs5UtA3AhlKqXSlVBAwE1hyyjlLgJ+oFmOAY64aPxdCCNG+cw65aK2tSqk5wGe0TFt8RWudo5S6vfX4AmA5LTNc8miZtvgz10U2nFcPGXWSvGff52/vF3zwPTs0D10IIYTnk7VchBDCR0hBF0IIHyEFvQuUUncrpbRSyqcXplBKPamU2qOU2q6U+kApFWN0JldRSk1SSuUqpfKUUvOMzuNqSqkUpdSXSqndSqkcpdRvjc7kLkqpAKXUltb7aHyCFPROUkql0LIcQqHRWdzgc2Cw1noILctA3GdwHpdoXeZiPjAZGATMUkoNMjaVy1mB32mtBwJjgF/7wXs+7rfAbqNDOJMU9M57Bvg97dxA5Wu01iu01tbWh9/Rcp+BLzqxzIXWugk4vsyFz9JalxxfSE9rXUNLgetlbCrXU0olA1cCLxmdxZmkoHeCUmoacEhrvc3oLAa4BfjE6BAu0gs42OZxEX5Q3I5TSqUBFwDrDY7iDs/S0iGzG5zDqWQ99DNQSn0B9Gzn0APA/cAV7k3kWmd7v1rrj1rPeYCWP9EXuTObGzm0hIUvUkpFAO8Dd2qtq43O40pKqanAUa31JqXUpQbHcSop6Gegtb68veeVUplAOrCtdVnOZGCzUmqU1vqwGyM61Zne73FKqZuBqcAE7bs3L/jlEhZKKTMtxXyR1vo/Rudxg7HAtNZlv0OAKKXUW1rrGw3O1WVyY1EXKaXygSyttTeu2uaQ1g1O/gZcorUuNTqPqyilAmn50HcCcIiWZS9u0FrnnPWFXqx1c5rXgQqt9Z0Gx3G71h763VrrqQZHcQoZQxeOeB6IBD5XSm1VSi0wOpArtH7we3yZi93A//lyMW81FrgJ+EHrf9utrT1X4YWkhy6EED5CeuhCCOEjpKALIYSPkIIuhBA+Qgq6EEL4CCnoQgjhI6SgCyGEj5CCLoQQPuL/AQ/FGXux1HZjAAAAAElFTkSuQmCC",
	"text/plain": [
	"<Figure size 432x288 with 1 Axes>"
	]
	},
	"metadata": {
	"needs_background": "light"
	},
	"output_type": "display_data"
	}
	],
	"source": [
	"bins = 50\n",
	"x_min = -5.0\n",
	"x_max = 5.0\n",
	"\n",
	"fig = plt.figure()\n",
	"ax = fig.add_subplot(111)\n",
	"ax.hist(samples, bins=bins, range=(x_min, x_max), density=True, color='tab:blue', label=\"samples\")\n",
	"xs = np.linspace(x_min, x_max, 1000)\n",
	"ys = inv_dist_fun(xs)\n",
	"ax.plot(xs, ys, color='tab:red', linestyle='--', label=\"target\")\n",
	"ax.set_ylim([0.0, 0.5])\n",
	"ax.legend(loc=\"upper right\")\n",
	"plt.show()"
	]
	},
	{
	"cell_type": "code",
	"execution_count": null,
	"metadata": {},
	"outputs": [],
	"source": []
	}
	],
	"metadata": {
	"interpreter": {
	"hash": "02d3c230780613a338b307efe1c8520401457e802ba0c889aabda605c4f627d4"
	},
	"kernelspec": {
	"display_name": "Python 3.8.11 64-bit ('base': conda)",
	"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.8.11"
	},
	"orig_nbformat": 4
	},
	"nbformat": 4,
	"nbformat_minor": 2
	}