tatsy/metropolis.ipynb

## metropolis.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 Metropolis(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",
    "\n",
    "    def sample(self):\n",
    "        x0 = self.x\n",
    "        x1 = self.rng.normal(self.x, self.sigma)\n",
    "        p0 = self.inv_dist_fun(x0)\n",
    "        p1 = self.inv_dist_fun(x1)\n",
    "        a = min(p1 / (p0 + 1.0e-12), 1.0)\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 = Metropolis(0.0, inv_dist_fun, 1.0)\n",
    "n_burnin = max(n_mutate // 100, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Burn-in\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": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuPUlEQVR4nO3deXyU1b3H8c+ZyWRfIQGBQBJ2kABCQBRELKUFZFGLC4p6ay2g17q02uJatfZet6q1opS61g2vooiIsisgCAn7jgFCEhIg+zbZZubcPxJogEAmycw8mZnf+/Xi9WLyPHnOd0Ly4+Q85zlHaa0RQgjh/UxGBxBCCOEaUtCFEMJHSEEXQggfIQVdCCF8hBR0IYTwEVLQhRDCRzhV0JVS45VSB5RS6UqpOY0cH6OUKlFKba//84TrowohhLiQgKZOUEqZgbnAOCAbSFVKLdZa7z3r1HVa60luyCiEEMIJzvTQhwPpWuvDWusaYAEw1b2xhBBCNFeTPXSgC5DV4HU2cGkj512mlNoB5AAPaq33nH2CUmomMBMgLCxsaN++fZufWAgh/NiWLVvytdZxjR1zpqCrRj529noBW4EErXW5UmoisAjodc4naT0fmA+QkpKi09LSnGheCCHEKUqpo+c75syQSzbQtcHreOp64adprUu11uX1f18KWJRSsS3IKoQQooWcKeipQC+lVJJSKhC4CVjc8ASl1EVKKVX/9+H11y1wdVghhBDn1+SQi9bappS6B1gGmIG3tdZ7lFKz64/PA6YBdymlbEAlcJOWZRyFEMKjlFF1V8bQhfAvtbW1ZGdnU1VVZXQUrxAcHEx8fDwWi+WMjyultmitUxr7HGduigohRKtlZ2cTERFBYmIi9SO04jy01hQUFJCdnU1SUpLTnyeP/gshPKKqqor27dtLMXeCUor27ds3+7cZKehCCI+RYu68lnytpKALIYSPkIIuhPALxcXFvP76625vZ9GiRezde/ZSV54hBV0I4ReaW9C11jgcjma3IwVdCCHcbM6cORw6dIjBgwfzwAMPMHbsWIYMGUJycjJffvklABkZGfTr14+7776bIUOGkJWVxV/+8hf69u3LuHHjmD59Oi+++CIAhw4dYvz48QwdOpQrrriC/fv3s2HDBhYvXsxDDz3E4MGDOXTokEffo0xbFEIY4uitt53zsYgJ42l38804KivJmjnrnONR115L9HXXYisq4ti9951xLOH9f1+wvWeffZbdu3ezfft2bDYbVquVyMhI8vPzGTFiBFOmTAHgwIEDvPPOO7z++uukpaWxcOFCtm3bhs1mY8iQIQwdOhSAmTNnMm/ePHr16sWmTZu4++67Wb16NVOmTGHSpElMmzatpV+aFpOCLoTwO1prHnnkEdauXYvJZOLYsWOcOHECgISEBEaMGAHA+vXrmTp1KiEhIQBMnjwZgPLycjZs2MD1119/+prV1dUefhfnkoIuhDDEhXrUppCQCx4PiIlpskd+IR9++CF5eXls2bIFi8VCYmLi6TnfYWFhp88735P0DoeD6Ohotm/f3uIM7iBj6EIIvxAREUFZWRkAJSUldOjQAYvFwpo1azh6tPEVaUeNGsVXX31FVVUV5eXlfP311wBERkaSlJTEp59+CtQV/h07dpzTjqdJQRdC+IX27dszcuRIBgwYwPbt20lLSyMlJYUPP/yQ8222M2zYMKZMmcKgQYO47rrrSElJISoqCqjr5b/11lsMGjSIiy+++PSN1ZtuuokXXniBSy65xOM3RWVxLiGER+zbt49+/foZHaPZysvLCQ8Px2q1Mnr0aObPn8+QIUM80nZjXzNZnEsIIVpo5syZ7N27l6qqKm6//XaPFfOWkIIuhBAX8NFHHxkdwWkyhi6EED5CCroQQvgIKehCCOEjpKALIYSPkJuiQghDJM752qXXy3j2apderzXGjBnDiy++SEpKo7ML3UZ66EII4SOkoAsh/EJFRQVXX301gwYNYsCAAXzyySc8/fTTDBs2jAEDBjBz5szTa7eMGTOGBx54gNGjR9OvXz9SU1O57rrr6NWrF4899hhQt9Ru3759uf322xk4cCDTpk3DarWe0+7y5cu57LLLGDJkCNdffz3l5eVA3XK+/fv3Z+DAgTz44IMueY9S0IUQfuHbb7+lc+fO7Nixg927dzN+/HjuueceUlNT2b17N5WVlSxZsuT0+YGBgaxdu5bZs2czdepU5s6dy+7du3n33XcpKCgA6pbanTlzJjt37iQyMvKcDTTy8/N55plnWLlyJVu3biUlJYWXXnqJwsJCvvjiC/bs2cPOnTtP/yfRWlLQhRB+ITk5mZUrV/KnP/2JdevWERUVxZo1a7j00ktJTk5m9erV7Nmz5/T5p9ZHT05O5uKLL6ZTp04EBQXRvXt3srKyAOjatSsjR44EYMaMGaxfv/6MNn/88Uf27t3LyJEjGTx4MO+99x5Hjx4lMjKS4OBg7rzzTj7//HNCQ0Nd8h7lpqgQwi/07t2bLVu2sHTpUh5++GF+8YtfMHfuXNLS0ujatStPPvnk6SV0AYKCggAwmUyn/37qtc1mA0ApdUYbZ7/WWjNu3Dg+/vjjc/Js3ryZVatWsWDBAl577TVWr17d6vcoPXQhhF/IyckhNDSUGTNm8OCDD7J161YAYmNjKS8v57PPPmv2NTMzM9m4cSMAH3/8MaNGjTrj+IgRI/jhhx9IT08HwGq1cvDgQcrLyykpKWHixIm88sorLltXXXroQghDeHqa4a5du3jooYcwmUxYLBbeeOMNFi1aRHJyMomJiQwbNqzZ1+zXrx/vvfces2bNolevXtx1111nHI+Li+Pdd99l+vTpp3c0euaZZ4iIiGDq1KlUVVWhtebll192yXuU5XOFEB7hrcvnnk9GRgaTJk1i9+7dbmujucvnypCLEEL4CCnoQgjRAomJiW7tnbeEFHQhhMcYNcTrjVrytZKCLoTwiODgYAoKCqSoO0FrTUFBAcHBwc36PJnlIoTwiPj4eLKzs8nLyzM6ilcIDg4mPj6+WZ8jBV0I4REWi4WkpCSjY/g0GXIRQggf4VRBV0qNV0odUEqlK6XmXOC8YUopu1JqmusiCiGEcEaTBV0pZQbmAhOA/sB0pVT/85z3HLDM1SGFEEI0zZke+nAgXWt9WGtdAywApjZy3u+AhcBJF+YTQgjhJGcKehcgq8Hr7PqPnaaU6gJcC8y70IWUUjOVUmlKqTS50y2EEK7lTEFXjXzs7ImkrwB/0lrbL3QhrfV8rXWK1jolLi7OyYjCSPbyCnKfeopKF60GJ4RwH2cKejbQtcHreCDnrHNSgAVKqQxgGvC6UuoaVwQUnle1bx+lK1YAUJuVSemXi8m4aTon//Y3tMNhcDohxPk4Mw89FeillEoCjgE3ATc3PEFrfXpyqVLqXWCJ1nqR62IKT6k5epTMO36DOSqKiDFjCO7Xj17r1nLiuecp+NebYDLT4YH7jY4phGhEkwVda21TSt1D3ewVM/C21nqPUmp2/fELjpsL76Frazn2wO9Ba7r+cx7KYgHAFBbGRU89CVpT8M9/EjJ4EBFXXWVsWCHEOZx6UlRrvRRYetbHGi3kWuv/an0sYYSCt9+hau9eurzyCoEJCWccU0px0eN1G9kGde9uRDwhRBPk0X8BgK2ggMxX57Kt0wD+8p0Nvvv6PGdeCm/sBnZ7fMcZIcSFSUEXAJijo5k76Dr2t0to8twOFYXcun8ZtsLLCGjXzgPphBDOkLVcBADKbGZVtxSOhTc9nTTIUcvPsrZS9MEHHkgmhHCWFHRB0aefkj//X+DkOtVZER1J7diXogWf4Kjf+FYIYTwp6H5O2+3kv/EGFT/8AKqxZ8gat6jnaOyFhZQuOd9YuxDC06Sg+7ny777DlpNLzM03N31yA9tjexLYowfFn33mpmRCiOaSgu7nihd+TkBcHBFjf9a8T1SKdrfeSsjAZLT9gis+CCE8RGa5+DF7SQnl69bR7pZbUAHN/1aIuelGN6QSQrSU9ND9mL2khPCRI4mcNKnF19AOB9bUVFnjRYg2QAq6Hwvs1o2u894gZMDFLb5G6TffcPTW26jatcuFyYQQLSEF3U/Zy8upPX681dcJHzUKAgIoW7nSBamEEK0hBd1PlS1bRvqYq6g+fLhV1zFHRRE2fDhlK1ainZzHLoRwDynofqps9RoCOnUiMCmp6ZObEDHu59RkZFBz6JALkgkhWkpmufghR3U1FRs2EH3tNahmPEx0tsQ5dQ8VtavUfAg88dDrfNJnrCzaJYRBpKD7oRtn/YOnKyu5+0g4W+e0/knPwpAo7rvyXg5HdXZBOiFES0lB90OXHt+LNSCIXbE9XHbNgzHdXHYtIUTLSEH3Qx/2HccPnZOpNbvunz+ktorpB1eyLa43IEMuQhhBCrofKgqOpCg40qXXrDZbGJ+xichqq0uvK4Rwnsxy8TPl69Yx6fAPmLRrn+x0mMxsi+vF0JP7ZfqiEAaRgu5nij75hF+lf4eDls9uOZ8tHfoQW1VK9U8/ufzaQoimSUH3I9pux7ppM9vjejVr7XNnbe3YB4CK9T+4/NpCiKZJQfcjVXv24CgrqyvobpAfEs2+mAR0TY1bri+EuDC5KepHKjb+CMD2uJ5ua+P3V/6OjNkyy0UII0gP3Y/U5uYQ1K8fJUERbm9LbowK4XlS0P1IpyefJOmTBW5tI8hWzaGJV1P47ntubUcIcS4p6H5GBQa69frVAUHgcGDdvNmt7QghziUF3U8U/vt9Mu/8Lbq21u1thQ4fjjUtTfYaFcLDpKD7ifL166g9nouyWNzeVujw4TjKyqjav9/tbQkh/kMKuh/QdjuV27YTOmSoR9oLHTYMAOsmGXYRwpOkoPuB6vR0HGVlhKZ4pqBbOnag3W/uIKhPb4+0J4SoI/PQ/YA1LQ2AEA/10AE6PvSQx9oSQtSRHrofCGjfnohf/hJLF89tQKG1pvrIEewlJR5rUwh/JwXdD0SOH0/8319p1XZzzVX9008cnjCRsjVrPNamEP5Ohlx81Kn9PoNsNSg0VQFBHm0/qGdPTOHhVG7bTvQ113i0bSH8lVM9dKXUeKXUAaVUulJqTiPHpyqldiqltiul0pRSo1wfVbTEyJxdfPb143Quz/Nou8pkImTwYCq3bfNou0L4syYLulLKDMwFJgD9gelKqf5nnbYKGKS1HgzcAbzp4pyihfoVZVBttnA8rL3H2w65ZDDVP/2EvazM420L4Y+cGXIZDqRrrQ8DKKUWAFOBvadO0FqXNzg/DJCVmdqI3kVZHIzuikN57nbJqeGewXm1/K/WXH/fv9jaoQ8Zz8oqjEK4kzM/5V2ArAavs+s/dgal1LVKqf3A19T10s+hlJpZPySTlpfn2SEAf2Sx19K9JIcDMd0MaX9fTDf+Mvx2DkZ3NaR9IfyNMwW9sakR5/TAtdZfaK37AtcAf2nsQlrr+VrrFK11SlxcXLOCiubrUZJDgHZwMMaYglodEMSGzsmUB4Ya0r4Q/saZgp4NNKwI8UDO+U7WWq8FeiilYluZTbRSXkg0/xwwhT3tEw3L0Kk8n2vTv3f5ptRCiHM5U9BTgV5KqSSlVCBwE7C44QlKqZ6qfpKzUmoIEAgUuDqsaJ6CkCgW9RztkQ0tzqdvUSYzd39FQulxwzII4S+avCmqtbYppe4BlgFm4G2t9R6l1Oz64/OAXwG3KaVqgUrgRi1b1hhu2PF9pEd3oSg40rAMe9slANC/MMOwDEL4C6ceLNJaLwWWnvWxeQ3+/hzwnGujidYIr7Hy9I9v8U7/Cfxf77GG5TgR2o7iwDB6F2U1fbIQolXk0X8f1bu4roAeNGiGy2lKcSCmG32KMo3NIYQfkILuo071iA9GxxucBA7GdKVzRT6Oykqjowjh02QtFx/VpyiTzPAOWC0hRkdhUY8r+KzXVRwMMT6LEL5Meug+SGtd94SoQfPPz2a1hFBjdv/Wd0L4O+mh+6gHR/83qg1NNLomfS15rx4m7t7fGR1FCJ8lPXQfpJQiNyyWnPC28zRu7+Isij/7zOgYQvg0Keg+qPSbb/hFxiajY5xhf0w3bCdPUnvihNFRhPBZUtB9UNFHHzPh6I9GxzjDqemTlTt3GpxECN8lBd3HaLudyj17OBCTYHSUMxyK6gwWC1VS0IVwGynoPqY6/RDaauVAG5nhckqt2ULY8OFGxxDCp8ksFx9TtauuB2z4E6KN6PaWbGQlhDtJD93H1GRnY4qKIseALeeEEMaSgu5jOtx/P72+W4P24JZzzrLl53N48mSKv1hkdBQhfFLb+6kXrWZqo4/Ym9u1o/b4CSq3bzc6ihA+SQq6D6nctZvMWbOoPnLE6CiNUiYTIckDqNwlM12EcAe5KerlEud8ffrv16SvZdbutYwIGQMGbmpxIcHJAyl4800clZVt9jcJIbyV9NB9SO/iLE6GRBu6Q1FTQgYNBLudqr17jY4ihM+Rgu5D+hZmtpkVFs8nZOBAIqdMxhQaanQUIXyODLn4iMjqCjpZC1iaNMLoKBcUEBtLl+efNzqGED5Jeug+IqKmgl3tu7OvXaLRUZqktaY2J8foGEL4HCnoPuJYRAf+eMXd7GmfZHSUJhW9/wHpPxuLrbDQ6ChC+BQp6D7CpB1GR3BacP9+gKy8KISrSUH3BVrz3rJnuHn/cqOTOCW4f38wmWTlRSFcTAq6D+hoLSS2qpTioHCjozjFFBpKUK9eVO7cZXQUIXyKFHQf0KcoC4ADbXCFxfMJGTiQyl270G1o31MhvJ1MW/QBvYsyqTYFkBHZyegoToue9ivCLr8M7HYIkG9DIVxBfpJ8QJ/iLNKj47GbzEZHuaCGyxSc9v0yMp692vNhhPBBUtB9wLrOg7AGBBkdo9l6FmcTYqs2OoYQPkMKug9Y3GOU0RFa5Le7FhPoqAXuNTqKED5Bbop6uThrEdFVZUbHaJEDMd3oUZKDo6bG6ChC+AQp6F7u5gMr+OeqF8ALZ4sciOmGxWGnev9+o6MI4ROkoHu5PkVZdSssKmV0lGY7Nc2ycoc8YCSEK0hB92KOigq6lR73qvnnDeWHRFEQHClLAAjhIlLQvVjl7j2Y0V5b0FGKRy//LRf9+c9GJxHCJ0hB92KVO3cA3vWE6NmORnbCHB5mdAwhfIJTBV0pNV4pdUApla6UmtPI8VuUUjvr/2xQSg1yfVRxtqiJE/mfYTMoDfLeghhRU0Heq69SuWu30VGE8HpNzkNXSpmBucA4IBtIVUot1lo33BTyCHCl1rpIKTUBmA9c6o7A4j8sXbqwrstgo2O0il2ZyH9jHgQEEJI8wOg4Qng1Z3row4F0rfVhrXUNsACY2vAErfUGrXVR/csfgXjXxhRnsxUVUfzZZ0RVe+cc9FOslhACu3enSlZeFKLVnCnoXYCsBq+z6z92Pr8BvmnsgFJqplIqTSmVlpeX53xKcQ5rWhq5jz1Opwrv3/UnZOBAKnfulJUXhWglZwp6YxOcG/3JU0pdRV1B/1Njx7XW87XWKVrrlLi4OOdTinNU7dwJFguHojobHaXVQgYmYy8spPaY7DMqRGs4U9Czga4NXscD5/zkKaUGAm8CU7XWBa6JJ86ncsdOgvv2pdZsMTpKqwUPHIgpNJTa7KymTxZCnJczBT0V6KWUSlJKBQI3AYsbnqCU6gZ8DtyqtT7o+piiIW23U7V7NyHJyUZHcYngfv3onbqZsBEjjI4ihFdrcpaL1tqmlLoHWAaYgbe11nuUUrPrj88DngDaA6+rukfQbVrrFPfF9m81R4/isFoJGTQQNhqdpvWUSR6HEMIVlFE3olJSUnRaWpohbfsCW2EhKjCIHs98Z3QUl7g8Zxc3HFzNQ1f8N7Xmun6GbHwhxLmUUlvO12GWrpGXCmjXzueesOxTnEXP4myjYwjhtWSDCy/RcPu2WTu/ZG/7BK9/qKihve0SAehfeIR97RMNzSKEt5IeupcJslUz+fB6EkuPGx3FpYqDIzgWFkv/ggyjowjhtaSge5lexdmY0ez34gW5zmdvu0QuLszwys06hGgLZMjFy/QtygTgQEyCwUlcL/WivgQ6bITYqqm0BBsdRwivIwXdy/QvyCA7PM6rV1g8n3VdBvvUfQEhPE0KupdRaHbE9jQ6hluF1lZitYQYHcMlGt7MbqjhlMya7GNUbPiB2sxMwseOJfSSSzwVT/gYKehe5qkRd/j0GPMftnxMn6JMZv680eWAfErN0aOceP4FylevBq1RFguhl9Y9LWvdupWCt96mw+8fIKhHD4OTCm8hN0W9kRduCO2srIgOdC3PI7K6wugobmUvr+DI9Tdg3bSJ2Ltm033pUvrs3EH4FaMAqD2Wg3XzZo5cex2FH3woK1EKp0gP3YvcteNzOlqLePKy3xgdxW1OzUfvV5hhaA536/HMd1ze91r2x3SjMCcK5u8F9p4eiomaPImwy0aQ+9jjnHjmGWqzMunwpz/JMgniguS7w4sMzk83OoLbHYzpRq0yc3HhEaOjuJ7WzNr5JVdmbwNgQ+dkCkOiznt6QGws8a/Ppd3tt1H43r8pWbz4vOcKAdJD9xrhNVa6lZ1kVdehRkdxqxqzhUPRXXzyAaNb9i/nmsPr+NQ8hu/jG7/x2fhN1IHsmvcG4aNHuzeg8HpS0L3EqSGIffVDEr7s015XYdIOphkdxIVG5O5mxoEVrOiawtv9m7/oWMSYMUDdjJjaY8cIu3S4ixMKXyAF3Uv0L8zApkwcjO7a9MlebkNn31jn/ZTO5Xk8uGUBB6Pj+cfgX7Xqpnbu449RtXcfSQsXEhh/oZ0ghT+Sgu4lMiI7sbj7KKoDAo2O4hGJJblYt24jdIj3z8m+LHcPNpOJvw6/rdU7THV66imO/GoaOX/4AwkffUjSo982ep4sPeyfpKB7ie/jLznvuKsvum/7p5x8fhWJCz42OkqrLew1huUJwygLbPnTvQ3H1q/qPZk/bvmIB258BHrKuLr4D5nl4gVsRUVE1Pj2vOyz7YjtSeXu3TgqvPd9Vx8+TK+iun1SW1PMz7Ym/hI2dezH7fu+oVN5vsuuK7yfFHQvUPzpZyxY+iThNVajo3jMzrgeYLNh3brN6Cgtoh0Och95lCd/fBuLvda1F1eKfwyexjeJIygLDHXttYVXk4LuBaypqWRFdKDcj35497RLBIsF6+ZNRkdpkZIvF1O5fTvv9J/Y6nHzxhSERDE/eapffU+IpklBb+O0zUblli3sjPWv9TyqA4IISU7GujnV6CjNZi8r4+Tf/kbIoEGs6ube5wZ6F2Xy+KZ3CXT1bwHCK0lBb+Oq9u7FYbWyK7a70VE8rtNfn6HrW28aHaPZ8ue+jr2ggI6PPYZW7v0RC7VVc3nubq5L/96t7QjvIAW9jbNu3gzArvb+1UMHCEpKwhwebnSMZjPHxBAzYwYhyQPc3tb2uF5s6DSA639aQ1R1udvbE22bTFts4yLGTyCgY0eKf/DP/3sL3n0XZTLT7rZbjY7itNhZMz3a3rv9J/DGqj3ccHA1/0qe4tG2Rdvin1XCiwTGdyFq8mSjYxjGuvFHCj/8wOgYTqn+6SfKVq70+FK3WREdWdltGJOObCDOWuTRtkXbIgW9DavJzKT4i0XYy/33V+mwUaOoPZpJzdGjRkdp0sm/vUTOI4/iMODf68O+43in/0SKgyI83rZoO6Sgt2FlK1aS+/DD6MpKo6MY5tSGD+Xr1xuc5MKsW7dR/t13tL/zTswRni+qeaExLOo5mlqzjKL6M/nXb8MqNm4ksHt3AuLijI5iGEtCApauXalYt552t9xidJxGaa3Je/llzLGxtJthbMarsrbQvSSHxDmNH5c1Xnyb9NDbKEd1NdbUVMJGjTQ6iqGUUkSMHQsB5ja7DVvFhg1YU1OJnT0bU6ixD/p0KzvBdelriS87aWgOYQzpobdRlVu3oqurCbv8cqOjGOY/C1INgPAB8PBSoO31MnVNDSEpQ4m+4Xqjo/BFjyuZemg9Nx1YyYspNxsdR3iYFPQ25lQRm/bTGm5TZoYtKaTq28Z2sfE/gfZaatzwGH1rRVx1FcnLrPDECqOjUBoUxpKky7ku/Xs+7juOY+H+O1znj2TIpY36rNdVzBj/BFUBQUZHaRPu2LOEeategDY07KLtdooXLsRRXW10lDN83vNKas0B3HhgldFRhIdJD70NKw1y3ZKr3i47vAOdrIX0KMkxOspppUuWkPvoY5giI42Ocobi4Aje6T+RvJAYo6MID5Meeht0Wc5uHtv0rt+tgX4hmy7qjx3F5bm7jI4C1I2b5/3jNYL796+7advGLO5xBRs7u3/pAdG2SEFvg0Yc383A/ENUWEKMjtJmlASFs7d9Epfl7jY6CgDFCxdSm51N3AP3o0xt88cotLaSW/Yto2NFgdFRhIc49Z2olBqvlDqglEpXSp0zw1Up1VcptVEpVa2UetD1Mf2H0g6GndjPlg59cLh5pT5vs6HTAJJKjxv61GjinK/p8+AX7HnuFXa1T2LAkpIztodrS0JsNdz402pu+GmN0VGEhzRZMZRSZmAuMAHoD0xXSvU/67RC4F7gRZcn9DO9i7KIqS5n80Vnf4nF+s4DeSP5GsxRUYbmiK4uIy8kivf6TQClDM1yIQUhUXybcCnjjqbKGi9+wpku4HAgXWt9WGtdAywApjY8QWt9UmudCsgq+600/Pg+7ChSO/Y1Okqbkx8azeIeozBHRxua40RYex4YfS97vGCN+k97XQXUTYMVvs+Zgt4FyGrwOrv+Y82mlJqplEpTSqXl5eW15BI+rzAkkhUJw2RrsfMIslVTvHAhNRkZhrQ/7PjeunXH23DPvKG80BhWdkthwtFNtK8sMTqOcDNnCnpj37ktmgystZ6vtU7RWqfE+fH6JBfyddLl/P2SG4yO0WYF22vIfeLPFC/83ONt1548ySOp73PHniUeb7s1Pun9M7Z06CPb1PkBZwp6NtC1wet4oO1MBvYhtcePY7HbjI7RppUERRA2aiQlS5agHQ6Ptp3/xhsEOOws6PNzj7bbWifC2vPUiDvIDY81OopwM2cKeirQSymVpJQKBG4CFrs3ln/KfeIJXv7+VaNjtHlRk6dgy8316AbSNUePUvzpZ3yTOILcMO8sjBdVFFD67TKjYwg3avJJUa21TSl1D7AMMANva633KKVm1x+fp5S6CEgDIgGHUup+oL/WutR90X2LvbSUio0/sj3Bfxfjctbw72p43xLC+4++wrPDZpz+uDsX7cr7+6soi4WPvax33tAt+5eTs34vocOHEdCundFxhBs4NdFZa71Ua91ba91Da/3X+o/N01rPq//7ca11vNY6UmsdXf93KebNULZqNdTWsrbLIKOjtHnVAYGs6DaM2MpiTNr9wy7aZkNrB+1//V8UBbetx/yb45PeY9FVVRS+867RUYSbyFoubUTpt99g6dyZg9Fdmz5Z8M7FE7GZPPPtqwICiH/55box+0e+8Uib7pAd0YHICRMo+vBD2t3xawJiZK0XXyOPIrYB9pISKjZsJGL8eK+ZDme0U8U8vMaKyWF3Wzvl69ZTnZ4O0GYf8W+O2Ltm47BaKfz3v42OItzA+79DfYApMpLED94n5ubpRkfxKkklOXzw7dOMPrbDLde3l5WRM2cOx598yi3XN0JQr15ETpqErq4xOopwAxlyaQOUUoQMOjV27p7i5IsyIi8iN6w9Nx5cxffxg11+/ZMvvYS9qIgO//yny69tpM4vPI+S3wR9kvTQDWYrKCD3iT9Tk5lpdBSvo5WJT3qPJbHsBJfl7nHptSfd+SrFHy9gYdIo+n2QQeKcr9vsIlzNdaqYW7dtw15WZnAa4UrSQzfIqeJwTfpaZu1ezI353ciMvMjgVN5nbZdBzNi/nBn7l6FtD6ICWv8t7aiq4r5tn5Ib2o73+/3SBSnbnpqjRzk6/Wbaz5pFhwfuNzqOcBHpoRvs55lpHIjuKsW8hRwmM+/0n0hC6QmsW7a65qJKsbHzAF4dPI1qH90CMDAhgcirr6bwvfeoPXnS6DjCRaSgG6h78TF6lOawoluK0VG82g+dk5k19iHCLh3e6mtprTEFBfH2xZPY3qG3C9K1XXH33Yu22cif+7rRUYSLSEE30LjMNGpNZr6Pv8ToKN5NKbIjOgBQ/dNPLb5MTWYmR665lspdbWNXJHcL7NaNmBtuoPizz6g+csToOMIFpKAbqMISzMquKbJUrouUr1vH4clTKP3222Z/rq2oiKyZs7AdP445ynufBm2u2LvvIqBDB2qkoPsEuSlqoA989IabUcJGjCB44EByH32MwKTuBPdxbsjEXlZG9uy7qD12jG7vvkNgt25A29iM2t0CYmPpuWK5S24mC+PJv6IBtNb0Lcxgf0yCPBnqQkmPLye20xRePvgqJ2+8jTmjZnMsPO6Ci3bZS0vJ/PUdVB08SJeX/kbo0KEeTOx555t6eeR/JlC2ciURY8eizGYPpxKuIkMuBqjcto2X177GmOxtRkfxOfkh0Tx+2Z0EOGw8v+51gm3VFzzfFBqKOSaG+H+8SuS4cR5K2faUr13LsXvvo/hzz28cIlxHCroBCt99jzJLCBs7XWx0FJ+UEdWZB0ffw4I+Y6mqn3ZYsXkztSdOUpuTQ+ny5WT/7l5seXmogAC6/ms+EWPGGBvaYOFXXknI0KHkvfQy9hLZqs5bSUH3sOrDhylbsYIlSZf77BzntuBYeBxfdR8FwMTfvkbmbbeTfuWVpP9sLMfuvY/MtRu56Y/vAchj8NR9DS567FHsJSXk/eM1o+OIFpIxdA8rmP8vVFAQX/a4wugofuNQVBcevnwmnSvy0SgyIzqyr10CDpPZZx7nd4Xgfv2IuelGij76iKgpkwkZONDoSKKZpKB7kKO6morNm4i+4XpKSsONjuM3qgMC2d6hN9vx7QeFXCHu97+n6sBBdI2sxuiNlNbakIZTUlJ0WlqaIW0byVFTg66upsdf1xodRQinuHNrP9F8SqktWutGHy+XHrqbnfqVPrqqjApLCLVm+ZKLts9ir+XmAytZ31mGXbyJ3BT1kPu3/R9/W/sPMOg3IiGaI8hey88zU5mT9gGOigqj4wgnSUH3gEtOHuTSE/tYGz9YHiQSXqE8MJTnU26hU3k+x5/+i9FxhJOkoLtZkK2G321fSHZYLF/WT6MTwhvsiu3BR33HUfLllxR/scjoOMIJUtDdbPqBlXSyFvCPwdOoNVuMjiNEsyzo83NChw/n5LPPYi+XoZe2Tu7QuZF2OOhbdJTl3VLYGdfT6DhCNJtDmejyysvU5uZiDg8zOo5oghR0N1ImE49cPpNAh83oKEK0WEC7dgS0awdAyZKvCR99BeZI/1li2JvIkIsbaK0peOttbIWFOEzm0+uJCOHNarKzyXn4YbJmzZaZL22UFHQ3yJ/7OidfeIHSJUuMjiKEywTGx9PlxRep3LGDrLv/W4p6GyQF3cWKFy0i/7XXiLrmGmJuvdXoOEK4VOQvf0HnZ/8Xa2oqR399B7aiIqMjiQZkDN1FEud8zVVZW/nDlo/ZFduDxx0jsD281OhYQrhc1JQpmMLCOPb7P1CxYQNRV8vSAG2FFHQXCXDYuOnASnbFdufJEXdgM8mXVviuiLFj6bHsWywXXQRAbW4ulk6dDE4lpOq0krbZ0HY7NlMAD4+cRYUlhOqAQKNjCeF2vV7ZAkBiSQ5///5VVnUdypsDJrH3b9MMTua/pKC3QvXhw+TMeZjAhAQwjaYwJMroSEK4XFNrxh8L78Di7iO5Nn0tKSf3U/pNGBHjx8vGIQaQm6ItYC8v5+Qrr3Dk2uuoPXqUiJ9dJWu0CL9Vaw7grQGT+cPoeyi3hHLsgd9z9NZb0Xa70dH8jvTQm2nyb/7Ow2kfEFNdzpr4S/jXgMkUfe8wOpYQhjvQLoF7rnqArZfasOXno8xmtNaULVtO2KhR8qSpB0hBb4KjupqK9esxR0cTOnQoueGxHIrqwvv9fsnBmG5GxxOiTXEoE9HXXXv6ddWuXRy7/35UYCDhV44m4he/JOyyEQTExhqY0nc5VdCVUuOBvwNm4E2t9bNnHVf1xycCVuC/tNZbXZzVYyo2bMC6ZSvWtDQqt29HV1cTOXEioUOHkh8SzeOX/9boiEJ4heDkZBI++ojSb74h/dMvab9iJQB/HHUXu2J7cFFFAav/62KCevfGHBFhcFrv1+QWdEopM3AQGAdkA6nAdK313gbnTAR+R11BvxT4u9b60gtd1xVb0GmtwW6vG6uz2zGFhgJ1Y9wOqxUcDrTNjq6uQtfWEty3LwAVGzdSk5GBLb8AW14etpMnMcfE0Pl//weAI9f9ioq9+zgc1Zndsd1J69iXHbE9sZvMrcorhD8zaQc9i7MZlJfOV91HUhUQxM37l3Pr/uV1xyMjsXTqhKVzZzq/8Dzm8HCsW7ZQk5GBKTISc0QEKigIU1AQwf37A2AvKwO7HQIsKJMCkwlMJkyBdTPNtNY+d3O2tVvQDQfStdaH6y+2AJgK7G1wzlTg37ruf4cflVLRSqlOWuvcVmZv1HMT7mTSkQ2Y+M9/RtWmAK6Z8iwZz17N8aefpnTxV2d8TnFgGNMnPgXA4z++w+XH9wBQFBROYXAkmREdeb7+bn7nzpMpTpqO1RLijvhC+CWHMnEwptsZQ5VLEy/jp+h4EspOEGctokNJMe2PH2T0X9Zw5LnJlHz1FcULPjnjOiowkL47dwBw4pm/UvLll2ccN8fE0HvjBgCy7/kd5atWnS70Siks3brR4+u6ZTkyZ83CumnzGZ8f1Ls3Sf9X1+bRGbdSuWfPGcdDBg8i4Z13gLrOX/WRI2ccDxt5OV1few2AQ+MnUHvy5BnHI8f9nM7PPefkV615nOmhTwPGa63vrH99K3Cp1vqeBucsAZ7VWq+vf70K+JPWOu2sa80EZta/7AMccNUb8aBYIN/oEB4m79n3+dv7Be99zwla67jGDjjTQ2/s95Wz/xdw5hy01vOB+U602WYppdLO9+uOr5L37Pv87f2Cb75nZ+ahZwNdG7yOB3JacI4QQgg3cqagpwK9lFJJSqlA4CZg8VnnLAZuU3VGACXuGj8XQgjRuCaHXLTWNqXUPcAy6qYtvq213qOUml1/fB6wlLoZLunUTVv8tfsiG86rh4xaSN6z7/O39ws++J6bvCkqhBDCO8haLkII4SOkoAshhI+Qgt4KSqkHlVJaKeXTC1MopV5QSu1XSu1USn2hlIo2OpO7KKXGK6UOKKXSlVJzjM7jbkqprkqpNUqpfUqpPUqp+4zO5ClKKbNSalv9czQ+QQp6CymlulK3HEKm0Vk8YAUwQGs9kLplIB42OI9b1C9zMReYAPQHpiul+hubyu1swB+01v2AEcB/+8F7PuU+YJ/RIVxJCnrLvQz8kUYeoPI1WuvlWmtb/csfqXvOwBedXuZCa10DnFrmwmdprXNPLaSntS6jrsB1MTaV+yml4oGrgTeNzuJKUtBbQCk1BTimtd5hdBYD3AF8Y3QIN+kCZDV4nY0fFLdTlFKJwCXAJoOjeMIr1HXIfGozA1kP/TyUUiuBixo59CjwCPALzyZyrwu9X631l/XnPErdr+gfejKbBzm1hIUvUkqFAwuB+7XWpUbncSel1CTgpNZ6i1JqjMFxXEoK+nlorX/e2MeVUslAErCjflnOeGCrUmq41vq4ByO61Pne7ylKqduBScBY7bsPL/jlEhZKKQt1xfxDrfXnRufxgJHAlPplv4OBSKXUB1rrGQbnajV5sKiVlFIZQIrW2htXbXNK/QYnLwFXaq3zjM7jLkqpAOpu+o4FjlG37MXNWus9F/xEL1a/Oc17QKHW+n6D43hcfQ/9Qa31JIOjuISMoQtnvAZEACuUUtuVUvOMDuQO9Td+Ty1zsQ/4P18u5vVGArcCP6v/t91e33MVXkh66EII4SOkhy6EED5CCroQQvgIKehCCOEjpKALIYSPkIIuhBA+Qgq6EEL4CCnoQgjhI/4fbJwfbEvtRN8AAAAASUVORK5CYII=",
      "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 Metropolis(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",
	"\n",
	" def sample(self):\n",
	" x0 = self.x\n",
	" x1 = self.rng.normal(self.x, self.sigma)\n",
	" p0 = self.inv_dist_fun(x0)\n",
	" p1 = self.inv_dist_fun(x1)\n",
	" a = min(p1 / (p0 + 1.0e-12), 1.0)\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 = Metropolis(0.0, inv_dist_fun, 1.0)\n",
	"n_burnin = max(n_mutate // 100, 100)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 6,
	"metadata": {},
	"outputs": [],
	"source": [
	"# Burn-in\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": 7,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuPUlEQVR4nO3deXyU1b3H8c+ZyWRfIQGBQBJ2kABCQBRELKUFZFGLC4p6ay2g17q02uJatfZet6q1opS61g2vooiIsisgCAn7jgFCEhIg+zbZZubcPxJogEAmycw8mZnf+/Xi9WLyPHnOd0Ly4+Q85zlHaa0RQgjh/UxGBxBCCOEaUtCFEMJHSEEXQggfIQVdCCF8hBR0IYTwEVLQhRDCRzhV0JVS45VSB5RS6UqpOY0cH6OUKlFKba//84TrowohhLiQgKZOUEqZgbnAOCAbSFVKLdZa7z3r1HVa60luyCiEEMIJzvTQhwPpWuvDWusaYAEw1b2xhBBCNFeTPXSgC5DV4HU2cGkj512mlNoB5AAPaq33nH2CUmomMBMgLCxsaN++fZufWAgh/NiWLVvytdZxjR1zpqCrRj529noBW4EErXW5UmoisAjodc4naT0fmA+QkpKi09LSnGheCCHEKUqpo+c75syQSzbQtcHreOp64adprUu11uX1f18KWJRSsS3IKoQQooWcKeipQC+lVJJSKhC4CVjc8ASl1EVKKVX/9+H11y1wdVghhBDn1+SQi9bappS6B1gGmIG3tdZ7lFKz64/PA6YBdymlbEAlcJOWZRyFEMKjlFF1V8bQhfAvtbW1ZGdnU1VVZXQUrxAcHEx8fDwWi+WMjyultmitUxr7HGduigohRKtlZ2cTERFBYmIi9SO04jy01hQUFJCdnU1SUpLTnyeP/gshPKKqqor27dtLMXeCUor27ds3+7cZKehCCI+RYu68lnytpKALIYSPkIIuhPALxcXFvP76625vZ9GiRezde/ZSV54hBV0I4ReaW9C11jgcjma3IwVdCCHcbM6cORw6dIjBgwfzwAMPMHbsWIYMGUJycjJffvklABkZGfTr14+7776bIUOGkJWVxV/+8hf69u3LuHHjmD59Oi+++CIAhw4dYvz48QwdOpQrrriC/fv3s2HDBhYvXsxDDz3E4MGDOXTokEffo0xbFEIY4uitt53zsYgJ42l38804KivJmjnrnONR115L9HXXYisq4ti9951xLOH9f1+wvWeffZbdu3ezfft2bDYbVquVyMhI8vPzGTFiBFOmTAHgwIEDvPPOO7z++uukpaWxcOFCtm3bhs1mY8iQIQwdOhSAmTNnMm/ePHr16sWmTZu4++67Wb16NVOmTGHSpElMmzatpV+aFpOCLoTwO1prHnnkEdauXYvJZOLYsWOcOHECgISEBEaMGAHA+vXrmTp1KiEhIQBMnjwZgPLycjZs2MD1119/+prV1dUefhfnkoIuhDDEhXrUppCQCx4PiIlpskd+IR9++CF5eXls2bIFi8VCYmLi6TnfYWFhp88735P0DoeD6Ohotm/f3uIM7iBj6EIIvxAREUFZWRkAJSUldOjQAYvFwpo1azh6tPEVaUeNGsVXX31FVVUV5eXlfP311wBERkaSlJTEp59+CtQV/h07dpzTjqdJQRdC+IX27dszcuRIBgwYwPbt20lLSyMlJYUPP/yQ8222M2zYMKZMmcKgQYO47rrrSElJISoqCqjr5b/11lsMGjSIiy+++PSN1ZtuuokXXniBSy65xOM3RWVxLiGER+zbt49+/foZHaPZysvLCQ8Px2q1Mnr0aObPn8+QIUM80nZjXzNZnEsIIVpo5syZ7N27l6qqKm6//XaPFfOWkIIuhBAX8NFHHxkdwWkyhi6EED5CCroQQvgIKehCCOEjpKALIYSPkJuiQghDJM752qXXy3j2apderzXGjBnDiy++SEpKo7ML3UZ66EII4SOkoAsh/EJFRQVXX301gwYNYsCAAXzyySc8/fTTDBs2jAEDBjBz5szTa7eMGTOGBx54gNGjR9OvXz9SU1O57rrr6NWrF4899hhQt9Ru3759uf322xk4cCDTpk3DarWe0+7y5cu57LLLGDJkCNdffz3l5eVA3XK+/fv3Z+DAgTz44IMueY9S0IUQfuHbb7+lc+fO7Nixg927dzN+/HjuueceUlNT2b17N5WVlSxZsuT0+YGBgaxdu5bZs2czdepU5s6dy+7du3n33XcpKCgA6pbanTlzJjt37iQyMvKcDTTy8/N55plnWLlyJVu3biUlJYWXXnqJwsJCvvjiC/bs2cPOnTtP/yfRWlLQhRB+ITk5mZUrV/KnP/2JdevWERUVxZo1a7j00ktJTk5m9erV7Nmz5/T5p9ZHT05O5uKLL6ZTp04EBQXRvXt3srKyAOjatSsjR44EYMaMGaxfv/6MNn/88Uf27t3LyJEjGTx4MO+99x5Hjx4lMjKS4OBg7rzzTj7//HNCQ0Nd8h7lpqgQwi/07t2bLVu2sHTpUh5++GF+8YtfMHfuXNLS0ujatStPPvnk6SV0AYKCggAwmUyn/37qtc1mA0ApdUYbZ7/WWjNu3Dg+/vjjc/Js3ryZVatWsWDBAl577TVWr17d6vcoPXQhhF/IyckhNDSUGTNm8OCDD7J161YAYmNjKS8v57PPPmv2NTMzM9m4cSMAH3/8MaNGjTrj+IgRI/jhhx9IT08HwGq1cvDgQcrLyykpKWHixIm88sorLltXXXroQghDeHqa4a5du3jooYcwmUxYLBbeeOMNFi1aRHJyMomJiQwbNqzZ1+zXrx/vvfces2bNolevXtx1111nHI+Li+Pdd99l+vTpp3c0euaZZ4iIiGDq1KlUVVWhtebll192yXuU5XOFEB7hrcvnnk9GRgaTJk1i9+7dbmujucvnypCLEEL4CCnoQgjRAomJiW7tnbeEFHQhhMcYNcTrjVrytZKCLoTwiODgYAoKCqSoO0FrTUFBAcHBwc36PJnlIoTwiPj4eLKzs8nLyzM6ilcIDg4mPj6+WZ8jBV0I4REWi4WkpCSjY/g0GXIRQggf4VRBV0qNV0odUEqlK6XmXOC8YUopu1JqmusiCiGEcEaTBV0pZQbmAhOA/sB0pVT/85z3HLDM1SGFEEI0zZke+nAgXWt9WGtdAywApjZy3u+AhcBJF+YTQgjhJGcKehcgq8Hr7PqPnaaU6gJcC8y70IWUUjOVUmlKqTS50y2EEK7lTEFXjXzs7ImkrwB/0lrbL3QhrfV8rXWK1jolLi7OyYjCSPbyCnKfeopKF60GJ4RwH2cKejbQtcHreCDnrHNSgAVKqQxgGvC6UuoaVwQUnle1bx+lK1YAUJuVSemXi8m4aTon//Y3tMNhcDohxPk4Mw89FeillEoCjgE3ATc3PEFrfXpyqVLqXWCJ1nqR62IKT6k5epTMO36DOSqKiDFjCO7Xj17r1nLiuecp+NebYDLT4YH7jY4phGhEkwVda21TSt1D3ewVM/C21nqPUmp2/fELjpsL76Frazn2wO9Ba7r+cx7KYgHAFBbGRU89CVpT8M9/EjJ4EBFXXWVsWCHEOZx6UlRrvRRYetbHGi3kWuv/an0sYYSCt9+hau9eurzyCoEJCWccU0px0eN1G9kGde9uRDwhRBPk0X8BgK2ggMxX57Kt0wD+8p0Nvvv6PGdeCm/sBnZ7fMcZIcSFSUEXAJijo5k76Dr2t0to8twOFYXcun8ZtsLLCGjXzgPphBDOkLVcBADKbGZVtxSOhTc9nTTIUcvPsrZS9MEHHkgmhHCWFHRB0aefkj//X+DkOtVZER1J7diXogWf4Kjf+FYIYTwp6H5O2+3kv/EGFT/8AKqxZ8gat6jnaOyFhZQuOd9YuxDC06Sg+7ny777DlpNLzM03N31yA9tjexLYowfFn33mpmRCiOaSgu7nihd+TkBcHBFjf9a8T1SKdrfeSsjAZLT9gis+CCE8RGa5+DF7SQnl69bR7pZbUAHN/1aIuelGN6QSQrSU9ND9mL2khPCRI4mcNKnF19AOB9bUVFnjRYg2QAq6Hwvs1o2u894gZMDFLb5G6TffcPTW26jatcuFyYQQLSEF3U/Zy8upPX681dcJHzUKAgIoW7nSBamEEK0hBd1PlS1bRvqYq6g+fLhV1zFHRRE2fDhlK1ainZzHLoRwDynofqps9RoCOnUiMCmp6ZObEDHu59RkZFBz6JALkgkhWkpmufghR3U1FRs2EH3tNahmPEx0tsQ5dQ8VtavUfAg88dDrfNJnrCzaJYRBpKD7oRtn/YOnKyu5+0g4W+e0/knPwpAo7rvyXg5HdXZBOiFES0lB90OXHt+LNSCIXbE9XHbNgzHdXHYtIUTLSEH3Qx/2HccPnZOpNbvunz+ktorpB1eyLa43IEMuQhhBCrofKgqOpCg40qXXrDZbGJ+xichqq0uvK4Rwnsxy8TPl69Yx6fAPmLRrn+x0mMxsi+vF0JP7ZfqiEAaRgu5nij75hF+lf4eDls9uOZ8tHfoQW1VK9U8/ufzaQoimSUH3I9pux7ppM9vjejVr7XNnbe3YB4CK9T+4/NpCiKZJQfcjVXv24CgrqyvobpAfEs2+mAR0TY1bri+EuDC5KepHKjb+CMD2uJ5ua+P3V/6OjNkyy0UII0gP3Y/U5uYQ1K8fJUERbm9LbowK4XlS0P1IpyefJOmTBW5tI8hWzaGJV1P47ntubUcIcS4p6H5GBQa69frVAUHgcGDdvNmt7QghziUF3U8U/vt9Mu/8Lbq21u1thQ4fjjUtTfYaFcLDpKD7ifL166g9nouyWNzeVujw4TjKyqjav9/tbQkh/kMKuh/QdjuV27YTOmSoR9oLHTYMAOsmGXYRwpOkoPuB6vR0HGVlhKZ4pqBbOnag3W/uIKhPb4+0J4SoI/PQ/YA1LQ2AEA/10AE6PvSQx9oSQtSRHrofCGjfnohf/hJLF89tQKG1pvrIEewlJR5rUwh/JwXdD0SOH0/8319p1XZzzVX9008cnjCRsjVrPNamEP5Ohlx81Kn9PoNsNSg0VQFBHm0/qGdPTOHhVG7bTvQ113i0bSH8lVM9dKXUeKXUAaVUulJqTiPHpyqldiqltiul0pRSo1wfVbTEyJxdfPb143Quz/Nou8pkImTwYCq3bfNou0L4syYLulLKDMwFJgD9gelKqf5nnbYKGKS1HgzcAbzp4pyihfoVZVBttnA8rL3H2w65ZDDVP/2EvazM420L4Y+cGXIZDqRrrQ8DKKUWAFOBvadO0FqXNzg/DJCVmdqI3kVZHIzuikN57nbJqeGewXm1/K/WXH/fv9jaoQ8Zz8oqjEK4kzM/5V2ArAavs+s/dgal1LVKqf3A19T10s+hlJpZPySTlpfn2SEAf2Sx19K9JIcDMd0MaX9fTDf+Mvx2DkZ3NaR9IfyNMwW9sakR5/TAtdZfaK37AtcAf2nsQlrr+VrrFK11SlxcXLOCiubrUZJDgHZwMMaYglodEMSGzsmUB4Ya0r4Q/saZgp4NNKwI8UDO+U7WWq8FeiilYluZTbRSXkg0/xwwhT3tEw3L0Kk8n2vTv3f5ptRCiHM5U9BTgV5KqSSlVCBwE7C44QlKqZ6qfpKzUmoIEAgUuDqsaJ6CkCgW9RztkQ0tzqdvUSYzd39FQulxwzII4S+avCmqtbYppe4BlgFm4G2t9R6l1Oz64/OAXwG3KaVqgUrgRi1b1hhu2PF9pEd3oSg40rAMe9slANC/MMOwDEL4C6ceLNJaLwWWnvWxeQ3+/hzwnGujidYIr7Hy9I9v8U7/Cfxf77GG5TgR2o7iwDB6F2U1fbIQolXk0X8f1bu4roAeNGiGy2lKcSCmG32KMo3NIYQfkILuo071iA9GxxucBA7GdKVzRT6Oykqjowjh02QtFx/VpyiTzPAOWC0hRkdhUY8r+KzXVRwMMT6LEL5Meug+SGtd94SoQfPPz2a1hFBjdv/Wd0L4O+mh+6gHR/83qg1NNLomfS15rx4m7t7fGR1FCJ8lPXQfpJQiNyyWnPC28zRu7+Isij/7zOgYQvg0Keg+qPSbb/hFxiajY5xhf0w3bCdPUnvihNFRhPBZUtB9UNFHHzPh6I9GxzjDqemTlTt3GpxECN8lBd3HaLudyj17OBCTYHSUMxyK6gwWC1VS0IVwGynoPqY6/RDaauVAG5nhckqt2ULY8OFGxxDCp8ksFx9TtauuB2z4E6KN6PaWbGQlhDtJD93H1GRnY4qKIseALeeEEMaSgu5jOtx/P72+W4P24JZzzrLl53N48mSKv1hkdBQhfFLb+6kXrWZqo4/Ym9u1o/b4CSq3bzc6ihA+SQq6D6nctZvMWbOoPnLE6CiNUiYTIckDqNwlM12EcAe5KerlEud8ffrv16SvZdbutYwIGQMGbmpxIcHJAyl4800clZVt9jcJIbyV9NB9SO/iLE6GRBu6Q1FTQgYNBLudqr17jY4ihM+Rgu5D+hZmtpkVFs8nZOBAIqdMxhQaanQUIXyODLn4iMjqCjpZC1iaNMLoKBcUEBtLl+efNzqGED5Jeug+IqKmgl3tu7OvXaLRUZqktaY2J8foGEL4HCnoPuJYRAf+eMXd7GmfZHSUJhW9/wHpPxuLrbDQ6ChC+BQp6D7CpB1GR3BacP9+gKy8KISrSUH3BVrz3rJnuHn/cqOTOCW4f38wmWTlRSFcTAq6D+hoLSS2qpTioHCjozjFFBpKUK9eVO7cZXQUIXyKFHQf0KcoC4ADbXCFxfMJGTiQyl270G1o31MhvJ1MW/QBvYsyqTYFkBHZyegoToue9ivCLr8M7HYIkG9DIVxBfpJ8QJ/iLNKj47GbzEZHuaCGyxSc9v0yMp692vNhhPBBUtB9wLrOg7AGBBkdo9l6FmcTYqs2OoYQPkMKug9Y3GOU0RFa5Le7FhPoqAXuNTqKED5Bbop6uThrEdFVZUbHaJEDMd3oUZKDo6bG6ChC+AQp6F7u5gMr+OeqF8ALZ4sciOmGxWGnev9+o6MI4ROkoHu5PkVZdSssKmV0lGY7Nc2ycoc8YCSEK0hB92KOigq6lR73qvnnDeWHRFEQHClLAAjhIlLQvVjl7j2Y0V5b0FGKRy//LRf9+c9GJxHCJ0hB92KVO3cA3vWE6NmORnbCHB5mdAwhfIJTBV0pNV4pdUApla6UmtPI8VuUUjvr/2xQSg1yfVRxtqiJE/mfYTMoDfLeghhRU0Heq69SuWu30VGE8HpNzkNXSpmBucA4IBtIVUot1lo33BTyCHCl1rpIKTUBmA9c6o7A4j8sXbqwrstgo2O0il2ZyH9jHgQEEJI8wOg4Qng1Z3row4F0rfVhrXUNsACY2vAErfUGrXVR/csfgXjXxhRnsxUVUfzZZ0RVe+cc9FOslhACu3enSlZeFKLVnCnoXYCsBq+z6z92Pr8BvmnsgFJqplIqTSmVlpeX53xKcQ5rWhq5jz1Opwrv3/UnZOBAKnfulJUXhWglZwp6YxOcG/3JU0pdRV1B/1Njx7XW87XWKVrrlLi4OOdTinNU7dwJFguHojobHaXVQgYmYy8spPaY7DMqRGs4U9Czga4NXscD5/zkKaUGAm8CU7XWBa6JJ86ncsdOgvv2pdZsMTpKqwUPHIgpNJTa7KymTxZCnJczBT0V6KWUSlJKBQI3AYsbnqCU6gZ8DtyqtT7o+piiIW23U7V7NyHJyUZHcYngfv3onbqZsBEjjI4ihFdrcpaL1tqmlLoHWAaYgbe11nuUUrPrj88DngDaA6+rukfQbVrrFPfF9m81R4/isFoJGTQQNhqdpvWUSR6HEMIVlFE3olJSUnRaWpohbfsCW2EhKjCIHs98Z3QUl7g8Zxc3HFzNQ1f8N7Xmun6GbHwhxLmUUlvO12GWrpGXCmjXzueesOxTnEXP4myjYwjhtWSDCy/RcPu2WTu/ZG/7BK9/qKihve0SAehfeIR97RMNzSKEt5IeupcJslUz+fB6EkuPGx3FpYqDIzgWFkv/ggyjowjhtaSge5lexdmY0ez34gW5zmdvu0QuLszwys06hGgLZMjFy/QtygTgQEyCwUlcL/WivgQ6bITYqqm0BBsdRwivIwXdy/QvyCA7PM6rV1g8n3VdBvvUfQEhPE0KupdRaHbE9jQ6hluF1lZitYQYHcMlGt7MbqjhlMya7GNUbPiB2sxMwseOJfSSSzwVT/gYKehe5qkRd/j0GPMftnxMn6JMZv680eWAfErN0aOceP4FylevBq1RFguhl9Y9LWvdupWCt96mw+8fIKhHD4OTCm8hN0W9kRduCO2srIgOdC3PI7K6wugobmUvr+DI9Tdg3bSJ2Ltm033pUvrs3EH4FaMAqD2Wg3XzZo5cex2FH3woK1EKp0gP3YvcteNzOlqLePKy3xgdxW1OzUfvV5hhaA536/HMd1ze91r2x3SjMCcK5u8F9p4eiomaPImwy0aQ+9jjnHjmGWqzMunwpz/JMgniguS7w4sMzk83OoLbHYzpRq0yc3HhEaOjuJ7WzNr5JVdmbwNgQ+dkCkOiznt6QGws8a/Ppd3tt1H43r8pWbz4vOcKAdJD9xrhNVa6lZ1kVdehRkdxqxqzhUPRXXzyAaNb9i/nmsPr+NQ8hu/jG7/x2fhN1IHsmvcG4aNHuzeg8HpS0L3EqSGIffVDEr7s015XYdIOphkdxIVG5O5mxoEVrOiawtv9m7/oWMSYMUDdjJjaY8cIu3S4ixMKXyAF3Uv0L8zApkwcjO7a9MlebkNn31jn/ZTO5Xk8uGUBB6Pj+cfgX7Xqpnbu449RtXcfSQsXEhh/oZ0ghT+Sgu4lMiI7sbj7KKoDAo2O4hGJJblYt24jdIj3z8m+LHcPNpOJvw6/rdU7THV66imO/GoaOX/4AwkffUjSo982ep4sPeyfpKB7ie/jLznvuKsvum/7p5x8fhWJCz42OkqrLew1huUJwygLbPnTvQ3H1q/qPZk/bvmIB258BHrKuLr4D5nl4gVsRUVE1Pj2vOyz7YjtSeXu3TgqvPd9Vx8+TK+iun1SW1PMz7Ym/hI2dezH7fu+oVN5vsuuK7yfFHQvUPzpZyxY+iThNVajo3jMzrgeYLNh3brN6Cgtoh0Och95lCd/fBuLvda1F1eKfwyexjeJIygLDHXttYVXk4LuBaypqWRFdKDcj35497RLBIsF6+ZNRkdpkZIvF1O5fTvv9J/Y6nHzxhSERDE/eapffU+IpklBb+O0zUblli3sjPWv9TyqA4IISU7GujnV6CjNZi8r4+Tf/kbIoEGs6ube5wZ6F2Xy+KZ3CXT1bwHCK0lBb+Oq9u7FYbWyK7a70VE8rtNfn6HrW28aHaPZ8ue+jr2ggI6PPYZW7v0RC7VVc3nubq5L/96t7QjvIAW9jbNu3gzArvb+1UMHCEpKwhwebnSMZjPHxBAzYwYhyQPc3tb2uF5s6DSA639aQ1R1udvbE22bTFts4yLGTyCgY0eKf/DP/3sL3n0XZTLT7rZbjY7itNhZMz3a3rv9J/DGqj3ccHA1/0qe4tG2Rdvin1XCiwTGdyFq8mSjYxjGuvFHCj/8wOgYTqn+6SfKVq70+FK3WREdWdltGJOObCDOWuTRtkXbIgW9DavJzKT4i0XYy/33V+mwUaOoPZpJzdGjRkdp0sm/vUTOI4/iMODf68O+43in/0SKgyI83rZoO6Sgt2FlK1aS+/DD6MpKo6MY5tSGD+Xr1xuc5MKsW7dR/t13tL/zTswRni+qeaExLOo5mlqzjKL6M/nXb8MqNm4ksHt3AuLijI5iGEtCApauXalYt552t9xidJxGaa3Je/llzLGxtJthbMarsrbQvSSHxDmNH5c1Xnyb9NDbKEd1NdbUVMJGjTQ6iqGUUkSMHQsB5ja7DVvFhg1YU1OJnT0bU6ixD/p0KzvBdelriS87aWgOYQzpobdRlVu3oqurCbv8cqOjGOY/C1INgPAB8PBSoO31MnVNDSEpQ4m+4Xqjo/BFjyuZemg9Nx1YyYspNxsdR3iYFPQ25lQRm/bTGm5TZoYtKaTq28Z2sfE/gfZaatzwGH1rRVx1FcnLrPDECqOjUBoUxpKky7ku/Xs+7juOY+H+O1znj2TIpY36rNdVzBj/BFUBQUZHaRPu2LOEeategDY07KLtdooXLsRRXW10lDN83vNKas0B3HhgldFRhIdJD70NKw1y3ZKr3i47vAOdrIX0KMkxOspppUuWkPvoY5giI42Ocobi4Aje6T+RvJAYo6MID5Meeht0Wc5uHtv0rt+tgX4hmy7qjx3F5bm7jI4C1I2b5/3jNYL796+7advGLO5xBRs7u3/pAdG2SEFvg0Yc383A/ENUWEKMjtJmlASFs7d9Epfl7jY6CgDFCxdSm51N3AP3o0xt88cotLaSW/Yto2NFgdFRhIc49Z2olBqvlDqglEpXSp0zw1Up1VcptVEpVa2UetD1Mf2H0g6GndjPlg59cLh5pT5vs6HTAJJKjxv61GjinK/p8+AX7HnuFXa1T2LAkpIztodrS0JsNdz402pu+GmN0VGEhzRZMZRSZmAuMAHoD0xXSvU/67RC4F7gRZcn9DO9i7KIqS5n80Vnf4nF+s4DeSP5GsxRUYbmiK4uIy8kivf6TQClDM1yIQUhUXybcCnjjqbKGi9+wpku4HAgXWt9WGtdAywApjY8QWt9UmudCsgq+600/Pg+7ChSO/Y1Okqbkx8azeIeozBHRxua40RYex4YfS97vGCN+k97XQXUTYMVvs+Zgt4FyGrwOrv+Y82mlJqplEpTSqXl5eW15BI+rzAkkhUJw2RrsfMIslVTvHAhNRkZhrQ/7PjeunXH23DPvKG80BhWdkthwtFNtK8sMTqOcDNnCnpj37ktmgystZ6vtU7RWqfE+fH6JBfyddLl/P2SG4yO0WYF22vIfeLPFC/83ONt1548ySOp73PHniUeb7s1Pun9M7Z06CPb1PkBZwp6NtC1wet4oO1MBvYhtcePY7HbjI7RppUERRA2aiQlS5agHQ6Ptp3/xhsEOOws6PNzj7bbWifC2vPUiDvIDY81OopwM2cKeirQSymVpJQKBG4CFrs3ln/KfeIJXv7+VaNjtHlRk6dgy8316AbSNUePUvzpZ3yTOILcMO8sjBdVFFD67TKjYwg3avJJUa21TSl1D7AMMANva633KKVm1x+fp5S6CEgDIgGHUup+oL/WutR90X2LvbSUio0/sj3Bfxfjctbw72p43xLC+4++wrPDZpz+uDsX7cr7+6soi4WPvax33tAt+5eTs34vocOHEdCundFxhBs4NdFZa71Ua91ba91Da/3X+o/N01rPq//7ca11vNY6UmsdXf93KebNULZqNdTWsrbLIKOjtHnVAYGs6DaM2MpiTNr9wy7aZkNrB+1//V8UBbetx/yb45PeY9FVVRS+867RUYSbyFoubUTpt99g6dyZg9Fdmz5Z8M7FE7GZPPPtqwICiH/55box+0e+8Uib7pAd0YHICRMo+vBD2t3xawJiZK0XXyOPIrYB9pISKjZsJGL8eK+ZDme0U8U8vMaKyWF3Wzvl69ZTnZ4O0GYf8W+O2Ltm47BaKfz3v42OItzA+79DfYApMpLED94n5ubpRkfxKkklOXzw7dOMPrbDLde3l5WRM2cOx598yi3XN0JQr15ETpqErq4xOopwAxlyaQOUUoQMOjV27p7i5IsyIi8iN6w9Nx5cxffxg11+/ZMvvYS9qIgO//yny69tpM4vPI+S3wR9kvTQDWYrKCD3iT9Tk5lpdBSvo5WJT3qPJbHsBJfl7nHptSfd+SrFHy9gYdIo+n2QQeKcr9vsIlzNdaqYW7dtw15WZnAa4UrSQzfIqeJwTfpaZu1ezI353ciMvMjgVN5nbZdBzNi/nBn7l6FtD6ICWv8t7aiq4r5tn5Ib2o73+/3SBSnbnpqjRzk6/Wbaz5pFhwfuNzqOcBHpoRvs55lpHIjuKsW8hRwmM+/0n0hC6QmsW7a65qJKsbHzAF4dPI1qH90CMDAhgcirr6bwvfeoPXnS6DjCRaSgG6h78TF6lOawoluK0VG82g+dk5k19iHCLh3e6mtprTEFBfH2xZPY3qG3C9K1XXH33Yu22cif+7rRUYSLSEE30LjMNGpNZr6Pv8ToKN5NKbIjOgBQ/dNPLb5MTWYmR665lspdbWNXJHcL7NaNmBtuoPizz6g+csToOMIFpKAbqMISzMquKbJUrouUr1vH4clTKP3222Z/rq2oiKyZs7AdP445ynufBm2u2LvvIqBDB2qkoPsEuSlqoA989IabUcJGjCB44EByH32MwKTuBPdxbsjEXlZG9uy7qD12jG7vvkNgt25A29iM2t0CYmPpuWK5S24mC+PJv6IBtNb0Lcxgf0yCPBnqQkmPLye20xRePvgqJ2+8jTmjZnMsPO6Ci3bZS0vJ/PUdVB08SJeX/kbo0KEeTOx555t6eeR/JlC2ciURY8eizGYPpxKuIkMuBqjcto2X177GmOxtRkfxOfkh0Tx+2Z0EOGw8v+51gm3VFzzfFBqKOSaG+H+8SuS4cR5K2faUr13LsXvvo/hzz28cIlxHCroBCt99jzJLCBs7XWx0FJ+UEdWZB0ffw4I+Y6mqn3ZYsXkztSdOUpuTQ+ny5WT/7l5seXmogAC6/ms+EWPGGBvaYOFXXknI0KHkvfQy9hLZqs5bSUH3sOrDhylbsYIlSZf77BzntuBYeBxfdR8FwMTfvkbmbbeTfuWVpP9sLMfuvY/MtRu56Y/vAchj8NR9DS567FHsJSXk/eM1o+OIFpIxdA8rmP8vVFAQX/a4wugofuNQVBcevnwmnSvy0SgyIzqyr10CDpPZZx7nd4Xgfv2IuelGij76iKgpkwkZONDoSKKZpKB7kKO6morNm4i+4XpKSsONjuM3qgMC2d6hN9vx7QeFXCHu97+n6sBBdI2sxuiNlNbakIZTUlJ0WlqaIW0byVFTg66upsdf1xodRQinuHNrP9F8SqktWutGHy+XHrqbnfqVPrqqjApLCLVm+ZKLts9ir+XmAytZ31mGXbyJ3BT1kPu3/R9/W/sPMOg3IiGaI8hey88zU5mT9gGOigqj4wgnSUH3gEtOHuTSE/tYGz9YHiQSXqE8MJTnU26hU3k+x5/+i9FxhJOkoLtZkK2G321fSHZYLF/WT6MTwhvsiu3BR33HUfLllxR/scjoOMIJUtDdbPqBlXSyFvCPwdOoNVuMjiNEsyzo83NChw/n5LPPYi+XoZe2Tu7QuZF2OOhbdJTl3VLYGdfT6DhCNJtDmejyysvU5uZiDg8zOo5oghR0N1ImE49cPpNAh83oKEK0WEC7dgS0awdAyZKvCR99BeZI/1li2JvIkIsbaK0peOttbIWFOEzm0+uJCOHNarKzyXn4YbJmzZaZL22UFHQ3yJ/7OidfeIHSJUuMjiKEywTGx9PlxRep3LGDrLv/W4p6GyQF3cWKFy0i/7XXiLrmGmJuvdXoOEK4VOQvf0HnZ/8Xa2oqR399B7aiIqMjiQZkDN1FEud8zVVZW/nDlo/ZFduDxx0jsD281OhYQrhc1JQpmMLCOPb7P1CxYQNRV8vSAG2FFHQXCXDYuOnASnbFdufJEXdgM8mXVviuiLFj6bHsWywXXQRAbW4ulk6dDE4lpOq0krbZ0HY7NlMAD4+cRYUlhOqAQKNjCeF2vV7ZAkBiSQ5///5VVnUdypsDJrH3b9MMTua/pKC3QvXhw+TMeZjAhAQwjaYwJMroSEK4XFNrxh8L78Di7iO5Nn0tKSf3U/pNGBHjx8vGIQaQm6ItYC8v5+Qrr3Dk2uuoPXqUiJ9dJWu0CL9Vaw7grQGT+cPoeyi3hHLsgd9z9NZb0Xa70dH8jvTQm2nyb/7Ow2kfEFNdzpr4S/jXgMkUfe8wOpYQhjvQLoF7rnqArZfasOXno8xmtNaULVtO2KhR8qSpB0hBb4KjupqK9esxR0cTOnQoueGxHIrqwvv9fsnBmG5GxxOiTXEoE9HXXXv6ddWuXRy7/35UYCDhV44m4he/JOyyEQTExhqY0nc5VdCVUuOBvwNm4E2t9bNnHVf1xycCVuC/tNZbXZzVYyo2bMC6ZSvWtDQqt29HV1cTOXEioUOHkh8SzeOX/9boiEJ4heDkZBI++ojSb74h/dMvab9iJQB/HHUXu2J7cFFFAav/62KCevfGHBFhcFrv1+QWdEopM3AQGAdkA6nAdK313gbnTAR+R11BvxT4u9b60gtd1xVb0GmtwW6vG6uz2zGFhgJ1Y9wOqxUcDrTNjq6uQtfWEty3LwAVGzdSk5GBLb8AW14etpMnMcfE0Pl//weAI9f9ioq9+zgc1Zndsd1J69iXHbE9sZvMrcorhD8zaQc9i7MZlJfOV91HUhUQxM37l3Pr/uV1xyMjsXTqhKVzZzq/8Dzm8HCsW7ZQk5GBKTISc0QEKigIU1AQwf37A2AvKwO7HQIsKJMCkwlMJkyBdTPNtNY+d3O2tVvQDQfStdaH6y+2AJgK7G1wzlTg37ruf4cflVLRSqlOWuvcVmZv1HMT7mTSkQ2Y+M9/RtWmAK6Z8iwZz17N8aefpnTxV2d8TnFgGNMnPgXA4z++w+XH9wBQFBROYXAkmREdeb7+bn7nzpMpTpqO1RLijvhC+CWHMnEwptsZQ5VLEy/jp+h4EspOEGctokNJMe2PH2T0X9Zw5LnJlHz1FcULPjnjOiowkL47dwBw4pm/UvLll2ccN8fE0HvjBgCy7/kd5atWnS70Siks3brR4+u6ZTkyZ83CumnzGZ8f1Ls3Sf9X1+bRGbdSuWfPGcdDBg8i4Z13gLrOX/WRI2ccDxt5OV1few2AQ+MnUHvy5BnHI8f9nM7PPefkV615nOmhTwPGa63vrH99K3Cp1vqeBucsAZ7VWq+vf70K+JPWOu2sa80EZta/7AMccNUb8aBYIN/oEB4m79n3+dv7Be99zwla67jGDjjTQ2/s95Wz/xdw5hy01vOB+U602WYppdLO9+uOr5L37Pv87f2Cb75nZ+ahZwNdG7yOB3JacI4QQgg3cqagpwK9lFJJSqlA4CZg8VnnLAZuU3VGACXuGj8XQgjRuCaHXLTWNqXUPcAy6qYtvq213qOUml1/fB6wlLoZLunUTVv8tfsiG86rh4xaSN6z7/O39ws++J6bvCkqhBDCO8haLkII4SOkoAshhI+Qgt4KSqkHlVJaKeXTC1MopV5QSu1XSu1USn2hlIo2OpO7KKXGK6UOKKXSlVJzjM7jbkqprkqpNUqpfUqpPUqp+4zO5ClKKbNSalv9czQ+QQp6CymlulK3HEKm0Vk8YAUwQGs9kLplIB42OI9b1C9zMReYAPQHpiul+hubyu1swB+01v2AEcB/+8F7PuU+YJ/RIVxJCnrLvQz8kUYeoPI1WuvlWmtb/csfqXvOwBedXuZCa10DnFrmwmdprXNPLaSntS6jrsB1MTaV+yml4oGrgTeNzuJKUtBbQCk1BTimtd5hdBYD3AF8Y3QIN+kCZDV4nY0fFLdTlFKJwCXAJoOjeMIr1HXIfGozA1kP/TyUUiuBixo59CjwCPALzyZyrwu9X631l/XnPErdr+gfejKbBzm1hIUvUkqFAwuB+7XWpUbncSel1CTgpNZ6i1JqjMFxXEoK+nlorX/e2MeVUslAErCjflnOeGCrUmq41vq4ByO61Pne7ylKqduBScBY7bsPL/jlEhZKKQt1xfxDrfXnRufxgJHAlPplv4OBSKXUB1rrGQbnajV5sKiVlFIZQIrW2htXbXNK/QYnLwFXaq3zjM7jLkqpAOpu+o4FjlG37MXNWus9F/xEL1a/Oc17QKHW+n6D43hcfQ/9Qa31JIOjuISMoQtnvAZEACuUUtuVUvOMDuQO9Td+Ty1zsQ/4P18u5vVGArcCP6v/t91e33MVXkh66EII4SOkhy6EED5CCroQQvgIKehCCOEjpKALIYSPkIIuhBA+Qgq6EEL4CCnoQgjhI/4fbJwfbEvtRN8AAAAASUVORK5CYII=",
	"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
	}