tatsy/simuanneal.ipynb

## simuanneal.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",
    "        self.t = 1.0\n",
    "        self.rho = 0.999\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) ** (1.0/self.t)\n",
    "        p1 = self.inv_dist_fun(x1) ** (1.0/self.t)\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",
    "        self.t *= self.rho\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/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtyklEQVR4nO3deXxU1f3/8ddnJpN9gSSELSzBQgFJQAiLQi2KWBAEtWrRitbWH6C11vUrVOvPWtufVqu2spW6YBV3XBBxA7WKKIR9R3YSFiEbZJksM3N+fyTwZQkwSWbmzvJ5Ph4+zMy9c897wuSTk3PPPVeMMSillAp9NqsDKKWU8g0t6EopFSa0oCulVJjQgq6UUmFCC7pSSoUJLehKKRUmvCroIjJCRLaIyDYRmdzA9qEiclhEVtf/95DvoyqllDqTqLPtICJ2YBowHCgA8kRknjFm40m7fm2MGe2HjEoppbzgTQ99ALDNGLPDGFMDvA6M9W8spZRSjXXWHjrQHsg/7nEBMLCB/c4XkTXAPuBeY8yGk3cQkQnABICEhIR+3bt3b3xiFXDr9h5u8Pns9ikBTqKUWrFiRaExplVD27wp6NLAcyevF7AS6GSMKReRy4D3gK6nvMiYWcAsgNzcXLN8+XIvmldW6zz5wwafX/7YqAAnUUqJyO7TbfNmyKUA6HDc40zqeuHHGGOOGGPK679eADhEJL0JWZVSSjWRNwU9D+gqIlkiEg2MA+Ydv4OItBERqf96QP1xi3wdViml1OmddcjFGOMSkduBTwA78IIxZoOITKrfPhO4GrhVRFyAExhndBlHpZQKKG/G0I8Ooyw46bmZx309FZjq22hKqXBSW1tLQUEBVVVVVkcJCbGxsWRmZuJwOLx+jVcFXSmlmqugoICkpCQ6d+5M/QitOg1jDEVFRRQUFJCVleX16/TSf6VUQFRVVZGWlqbF3AsiQlpaWqP/mtGCrpQKGC3m3mvK90oLulJKhQkt6EqpiFBaWsr06dP93s57773Hxo0nL3UVGFrQlVIRobEF3RiDx+NpdDta0JVSys8mT57M9u3b6dOnD3fddRfDhg2jb9++ZGdn8/777wOwa9cuevTowW233Ubfvn3Jz8/nz3/+M927d2f48OFcd911PPnkkwBs376dESNG0K9fP37yk5+wefNmlixZwrx587jvvvvo06cP27dvD+h71GmLSilL7B5/4ynPJY0cQer11+NxOsmfMPGU7SlXXkmLq67EVVLC3jt+f8K2Ti//54ztPfbYY6xfv57Vq1fjcrmorKwkOTmZwsJCBg0axJgxYwDYsmULL774ItOnT2f58uXMnTuXVatW4XK56Nu3L/369QNgwoQJzJw5k65du7J06VJuu+02Pv/8c8aMGcPo0aO5+uqrm/qtaTIt6EqpiGOM4Q9/+ANfffUVNpuNvXv38sMPPwDQqVMnBg0aBMDixYsZO3YscXFxAFx++eUAlJeXs2TJEq655ppjx6yurg7wuziVFnSllCXO1KO2xcWdcXtUy5Zn7ZGfyZw5czh06BArVqzA4XDQuXPnY3O+ExISju13uhVMPB4PLVq0YPXq1U3O4A86hq6UighJSUmUlZUBcPjwYTIyMnA4HHzxxRfs3t3wirRDhgzhgw8+oKqqivLycj78sG4p6eTkZLKysnjrrbeAusK/Zs2aU9oJNC3oSqmIkJaWxuDBg+nVqxerV69m+fLl5ObmMmfOHE53s53+/fszZswYevfuzVVXXUVubi4pKXU3dpkzZw7PP/88vXv35txzzz12YnXcuHE88cQTnHfeeQE/KSpWLYqoN7gIHae7wcUuvcGFaoRNmzbRo0cPq2M0Wnl5OYmJiVRWVnLhhRcya9Ys+vbtG5C2G/qeicgKY0xuQ/vrGLpSSp3BhAkT2LhxI1VVVdx0000BK+ZNoQVdKaXO4NVXX7U6gtd0DF0ppcKEFnSllAoTWtCVUipMaEFXSqkwoSdFlVKWON102KYKpmm0Q4cO5cknnyQ3t8HZhX6jPXSllAoTWtCVUhGhoqKCUaNG0bt3b3r16sUbb7zBI488Qv/+/enVqxcTJkw4tnbL0KFDueuuu7jwwgvp0aMHeXl5XHXVVXTt2pUHH3wQqFtqt3v37tx0003k5ORw9dVXU1lZeUq7n376Keeffz59+/blmmuuoby8HKhbzrdnz57k5ORw7733+uQ9akFXSkWEjz/+mHbt2rFmzRrWr1/PiBEjuP3228nLy2P9+vU4nU7mz59/bP/o6Gi++uorJk2axNixY5k2bRrr169n9uzZFBUVAXVL7U6YMIG1a9eSnJx8yg00CgsLefTRR1m4cCErV64kNzeXp556iuLiYt599102bNjA2rVrj/2SaC4t6EqpiJCdnc3ChQu5//77+frrr0lJSeGLL75g4MCBZGdn8/nnn7Nhw4Zj+x9dHz07O5tzzz2Xtm3bEhMTQ5cuXcjPzwegQ4cODB48GIAbbriBxYsXn9Dmd999x8aNGxk8eDB9+vThpZdeYvfu3SQnJxMbG8stt9zCO++8Q3x8vE/eo54UVcf4+iSVUsGkW7durFixggULFjBlyhQuvfRSpk2bxvLly+nQoQMPP/zwsSV0AWJiYgCw2WzHvj762OVyASAiJ7Rx8mNjDMOHD+e11147Jc+yZctYtGgRr7/+OlOnTuXzzz9v9nvUHrpSKiLs27eP+Ph4brjhBu69915WrlwJQHp6OuXl5bz99tuNPuaePXv49ttvAXjttdcYMmTICdsHDRrEN998w7Zt2wCorKzk+++/p7y8nMOHD3PZZZfxzDPP+Gxdde2hK6UsEehphuvWreO+++7DZrPhcDiYMWMG7733HtnZ2XTu3Jn+/fs3+pg9evTgpZdeYuLEiXTt2pVbb731hO2tWrVi9uzZXHfddcfuaPToo4+SlJTE2LFjqaqqwhjD008/7ZP3qMvnqmMaO+QSTPN+VfAL1eVzT2fXrl2MHj2a9evX+62Nxi6fq0MuSikVJrSgK6VUE3Tu3NmvvfOm0IKulAoYq4Z4Q1FTvlda0JVSAREbG0tRUZEWdS8YYygqKiI2NrZRr9NZLkqpgMjMzKSgoIBDhw5ZHSUkxMbGkpmZ2ajXaEFXSgWEw+EgKyvL6hhhTYdclFIqTHhV0EVkhIhsEZFtIjL5DPv1FxG3iFztu4hKKaW8cdaCLiJ2YBowEugJXCciPU+z3+PAJ74OqZRS6uy86aEPALYZY3YYY2qA14GxDez3O2AucNCH+ZRSSnnJm4LeHsg/7nFB/XPHiEh74Epg5pkOJCITRGS5iCzXM91KKeVb3hR0aeC5kyeSPgPcb4xxn+lAxphZxphcY0xuq1atvIyorBRXW8Vv18yle/Fuq6Mopc7Cm2mLBUCH4x5nAvtO2icXeL1+LeB04DIRcRlj3vNFSBVYXUr30qaymCXtsmlbUcTF+SsZvfNb3ux6EbN7jsSITo5SKhh5U9DzgK4ikgXsBcYB1x+/gzHm2ORSEZkNzNdiHpralhfy1yWzKIuOZ2mbnuxo0Z5fjniI/7P+A67d+gUesfFSz5FWx1RKNeCsBd0Y4xKR26mbvWIHXjDGbBCRSfXbzzhurkKH3eNmSt7LgOGhQb/BbbMDUBUVw7O9f44Yw7jvF7EptRPL2pwy0UkpZTGvrhQ1xiwAFpz0XIOF3Bjzq+bHUlb4+bYv6Xp4L3/pP579ieknbhRheu8rAchPzLAgnVLqbPTSfwWAq6iIcVsWsaRtLxa3793wPrYo/nneNQFOppTylhZ0BYC9RQum9b6KzamdzrpvRkUx4zd/gqv4fKJSUwOQTinlDZ2uoAAQu51FHXPZm3j26aQxnlouzl9JySuvBCCZUspbWtAVJW+9ReGsf4OX61TnJ7Umr3V3Sl5/A0/9jW+VUtbTgh7hjNtN4YwZVHzzDUhD15A17L0fXYi7uJgj8xt3Y2mllP9oQY9w5V9+iWvfflpef/3Zdz7O6vQfEX3OOZS+/bafkimlGksLeoQrnfsOUa1akTTs4sa9UITU8eOJy8nGuM+44oNSKkB0lksEcx8+TPnXX5P6y18iUY3/KLQc9ws/pFJKNZX20COY+/BhEgcPJnn06CYfw3g8VOblYTweHyZTSjWFFvQIFt2xIx1mziCu17lNPsaRjz5i9/gbqVq3zofJlFJNoQU9QrnLy6k9cKDZx0kcMgSioihbuNAHqZRSzaEFPUKVffIJ24ZeRPWOHc06jj0lhYQBAyj7bCHGy3nsSin/0IIeoco+/4Kotm2Jzso6+85nkTT8Emp27aJm+3YfJFNKNZUW9Ajkqa6mYskSki4aijTiYqLTSbx4GABlCxc1+1hKqabTaYsRqPK77zBOJ4kXNXLu+Wk4WmfQ+a03if3xj31yPKVU02hBj0BlX3yBLT6e+IEDfHbMuOxsnx1LKdU0OuQSgdJvu432z/4TW3S0z47pLq/g4JNPUv7NNz47plKqcbSHHoEcGRk4Mnx71yFbbAwlb76Fq6SExMGDfXpspZR3tIceYcq//priOXN8vv6KREWRcMEFVHy9WKcvKmURLegRpuSNNyh+4UWw+f6fPnHIYFwHD1K9davPj62UOjst6BHEuN1ULl1G/PmDfDJd8WQJQ4YAULFYx9GVsoIW9AhStWEDnrIyEs4/3y/Hd7RpQ1zv3piaGr8cXyl1ZnpSNIJUfPsdAAmDBvmtjc5vvO63Yyulzkx76BGkdv8+Ynr0ICotze9t6YlRpQJPC3oEafvww2T5uQftqaxk+2WjKJ79kl/bUUqdSgt6hBEfXkzUEFt8PHg8VC5b5td2lFKn0oIeIYr/8zJ7bvk/mNpav7cVP2AAlcuX671GlQowLegRonzx19Qe2I84HH5vK37AADxlZVRt3uz3tpRS/0sLegQwbjfOVauJ79svIO3F9+8PQOVSHXZRKpC0oEeA6m3b8JSVEZ8bmILuaJ1B6m9+TcyPuwWkPaVUHZ2HHgEqly8HIC5APXSA1vfdF7C2lFJ1tIceAaLS0kj62c9wtG8XsDaNMVTv3In78OGAtalUpNOCHgGSR4wg8x/P+GX9ltOp3rqVHSMvo+yLLwLWplKRTgt6mPM4nXgqKgLebsyPfoQtMRHnqtUBb1upSOVVQReRESKyRUS2icjkBraPFZG1IrJaRJaLyBDfR1VNUfbpp2wZMJCaXbsC2q7YbMT16YNz1aqAtqtUJDtrQRcROzANGAn0BK4TkZ4n7bYI6G2M6QP8GnjOxzlVE1WuXo0tNhZHhw4BbzvuvD5Ub92Ku6ws4G0rFYm86aEPALYZY3YYY2qA14Gxx+9gjCk3/7saUwKgKzMFiaq164jNzkbs9oC3Hd+3LxiDc83agLetVCTypqC3B/KPe1xQ/9wJRORKEdkMfEhdL/0UIjKhfkhm+aFDh5qSVzWCp7qaqi1biMvOtqT9uN69af/sP4nL7mVJ+0pFGm8KekNTI07pgRtj3jXGdAeuAP7c0IGMMbOMMbnGmNxWrVo1KqhqvOpNm8DlIjbHmoJui48nefhw7CkplrSvVKTxpqAXAMcPwGYC+063szHmK+AcEUlvZjbVTFFt29J6yuS6oQ+L1OzeTdGLs3WhLqUCwJuCngd0FZEsEYkGxgHzjt9BRH4k9ZOcRaQvEA0U+TqsahxH69ak3nRTQG5ocTrOtWs5+PjjeuNopQLgrAXdGOMCbgc+ATYBbxpjNojIJBGZVL/bz4H1IrKauhkxvzB6yxrLlf/3v7gsPlcRd955ADp9UakA8GoeujFmgTGmmzHmHGPMX+qfm2mMmVn/9ePGmHONMX2MMecbYxb7M7Q6O3dpKfkTJ1H6zruW5nC0b489NRXnuvWW5lAqEuiVomHKuX4DAHEWnRA9SkSIy87GuXaNpTmUigRa0MNU1bq6ud+xvayfMhibk03t7j14nE6roygV1rSghynn2nVEd+mCPSnJ6iik3ngj3fKWYYuLszqKUmFNC3oYMsbgXLfOsguKTmZPSsIWG2t1DKXCnt7gIkx1fnUOeDxWxzimaPZsPEfKaHXH76yOolTY0h56GBIRojt2JLpzZ6ujHFO1fgOlb79tdQylwpoW9DB05KOPgq54xuXk4Dp4kNoffrA6ilJhSwt6GCp59TVK3nzL6hgnODp90rlWV15Uyl+0oIcZ43bj3LCBuJwcq6OcIKZHD3A4qNKCrpTfaEEPM9XbtmMqKy2/oOhktpgYEgYMsDqGUmFNZ7mEmWMXFAXJlMXjdXxeb2SllD9pDz3M1BQUYEtJIbpTJ6ujKKUCTAt6mMm48066fvkFYgu+f1pXYSE7Lr+c0nffszqKUmEp+H7qVbMF6yX29tRUag/8gHP1aqujKBWWtKCHEee69eyZOJHqnTutjtIgsdmIy+6Fc53OdFHKH7SghxHnyhVU/PcrbAkJVkc5rdjsHKq3fK8rLyrlB1rQw4hz7Tqi2rbFkZFhdZTTiuudA243VRs3Wh1FqbCjBT2MONeuDZoVFk8nLieH5DGXY4uPtzqKUmFHC3qYcJWUUJufH3QXFJ0sKj2d9n/7G7E9elgdRamwowU9TLhLSonPzT12U+ZgZoyhdt8+q2MoFXa0oIeJmC5ZdHrlZeL79bM6ylmVvPwK2y4ehqu42OooSoUVLehhwrjdVkfwWmzPuuEWXXlRKd/Sgh4GjDFsu+hiDk2dZnUUr8T27Ak2m668qJSPaUEPA7V79+I6eJCo9DSro3jFFh9PTNeuONeuszqKUmFFC3oYONrTDcYVFk8nLicH57p1GGOsjqJU2NDlc8OAc+06JCaG2G7drI7itRZX/5yEC84Htxui9GOolC/oT1IYcK5bR2zPnojDYXUUr8X17k1c795Wx1AqrGhBDwPJI0YE9fotp+NcvwFPZYXeyUgpH9GCHgZSx99gdYQmOfj443iqq8l68w2roygVFvSkaIir3bcPV2Gh1TGaJDYnm+pNm/DU1FgdRamwoAU9xB2aPp0do0aH5GyRuJzemNpaqjdvtjqKUmFBC3qIq1q7jtjsbETE6iiNdnQhMecavcBIKV/QMfQQ5qmooHrbNpIuucSS9jtP/rDB53c9Nsqr10e1aUNUq1a6BIBSPqIFPYQ5128Aj6fuphEhSETo8PxzONq1tzqKUmFBC3oIc65dA0BsTmgWdCCkLoZSKth5VdBFZATwD8AOPGeMeeyk7b8E7q9/WA7caoxZ48ug6lQpl11GdGYmUS1bNup1pxsqsYKrpISSl18m8aKLicvuZXUcpULaWQu6iNiBacBwoADIE5F5xpjjbwq5E/ipMaZEREYCs4CB/gis/pejfXsc7UN7uEKioiicMROiorSgK9VM3sxyGQBsM8bsMMbUAK8DY4/fwRizxBhTUv/wOyDTtzHVyVwlJZS+/TauoiKrozSLPSmJ6C5dqNKVF5VqNm8Kensg/7jHBfXPnc5vgI8a2iAiE0RkuYgsP3TokPcp1Skqly9n/4N/pDY//+w7B7m4nByca9eG5Fx6pYKJNwW9oQnODf7kichF1BX0+xvaboyZZYzJNcbktmrVyvuU6hRVa9eCw0FMGNxsOS4nG3dxMbV79T6jSjWHNwW9AOhw3ONM4JSfPBHJAZ4DxhpjQnscIAQ416wltnt3bDExVkdptticHGzx8dQWhP5fG0pZyZtZLnlAVxHJAvYC44Drj99BRDoC7wDjjTHf+zylOoFxu6lav56UK66wOopPxPboQbe8ZYjdbnUUpULaWQu6McYlIrcDn1A3bfEFY8wGEZlUv30m8BCQBkyvvwTdZYzJ9V/syFazezeeysqQvaDoZGLTFSiU8gWv5qEbYxYAC056buZxX98C3OLbaOp0Yrp0oeuSb5Do0B9uOerIZ59RNOvfdJrzCrboaKvjKBWStGsUoqJSU7Enht5NLc6kat06qtZvsDqGUiFLC3oIOvDXv3LkowZnhoas+PPOA8C5aqXFSZQKXVrQQ4ynspKSV+ZQvXWr1VF8Kio9HUenjlSuXGV1FKVClhb0EFO1oW6FxVBekOt04s/ri3PlSr3ASKkm0tUWQ4xzTd2aZ3G9e1ucxPcSf3ohpqYaT0Vl2J0fUCoQtKCHmMqVq4ju3LnRKyyGguSRI0keOdLqGEqFLC3oocYY4geF90KW7rIy7ElJVscImJqCvVQs+YbaPXtIHDbs2AlipRpLC3qI6TBjeliPMe+7fzLOtWs556MFZ985xNXs3s0Pf3uC8s8/B2MQh4P4gYMAqFy5kqLnXyDj7ruIOecci5OqUKEnRUNQKN4Q2lvRXbpQs3MnrpKSs+8cwtzlFey85loqly4l/dZJdFmwgB+vXUPiT4YAULt3H5XLlrHzyqsofmVOWP8SV76jBT2EHHjkz+RPutXqGH4V3/fofPTwnr5oT0yg7aN/psuHH9LqjjuI6ZJ1wi/qlMtHc85HC0i44AJ+ePRRDj72GMbjsTCxCgVa0ENIxdKlEOY9tdjsbHA4cK4MvwuMjDEc+OtfOTy/7haAyZdeiqN1xmn3j0pPJ3P6NFJvupHil/7D4XnzAhVVhSgt6CHCXVpKzfbtxPXta3UUv7LFxhLbs0dYXmBUOHUaJf95merNm7x+jdhstJ4yhcyZM0gZM8aP6VQ40JOiIaJy9WoA4s7rY2mOQEj7zW/AE15/iZQtWkThtGmkXHEFre65p9GvTxo6FKibEVO7dy8JAwf4OKEKB1rQQ4Rz5ar6GylnWx3F75IvvdTqCD5Vs2sX++6fTOy559LmTw8366T2/j8+SNXGTWTNnUt0ZmjfIFz5ng65hIiYbt1I/eUvscXFWR0lIKq2fB82wy5liz5HoqLI/Oc/mn2HqbZ/+hN4POy75x6M2+2jhCpcaEEPESmjR9F6ymSrYwTM/of+yMG//c3qGD6R9ptf0+WjBTjaN79HHd2xI20e+iPONWsofvllH6RT4UQLeghwlZSE/bzskyUMGIhz/Xo8FRVWR2my6h07cK5bD+DTpRqSR48mcehQDj3zD2p27/bZcVXo04IeAkrfeputFwzGXVpqdZSAiR84EFyukB12MR4P+//wAPm33YqnutqnxxYR2vzpYVr+4lrsKSk+PbYKbVrQQ0BlXh7R53TB3qKF1VECJr7veeBwULlsqdVRmuTw+/Nwrl5Nxl13N3vcvCGO1q1pPWVKRH0m1NlpQQ9yxuXCuWIFCQMia5qaLT6euOxsKpflWR2l0dxlZRz8+9+J692blCvG+rUt57p15N9+O56qKr+2o0KDFvQgV7VxI57KSuL797c6SsC1/cujdHj+OatjNFrhtOm4i4po/eCDiM2/P2Ke8nLKFy6i+MUX/dqOCg1a0INc5bJlABFZ0GOysrAnJlodo9HsLVvS8oYbiMvu5fe2Es4/n8RLhlH03PO4iov93p4KbnphUZBLGjGSqNatiUpPtzqKJYpmz0ZsdlJvHG91FK+lT5wQ0PYy7rqLHZ+PoehfsyJqaqs6lfbQg1x0ZntSLr/c6hiWqfz2O4rnvGJ1DK9Ub91K2cKFAV/qNuacc0i58gpKXnuN2n37Atq2Ci5a0INYzZ49lL77Hu7ycqujWCZhyBBqd+8JifnWB//+FPv+8AAeC/69Wv32t2Tcczf2CP1LTtXRgh7Eyj5byP4pUzBOp9VRLHP0hg/lixdbnOTMKleuovzLL0m75RZLbp/naNeO1JtuwhYdHfC2VfDQgh7EKr79luguXYhq1crqKJZxdOqEo0MHKr4O3oJujOHQ009jT08n9YZfWprl8Lx5/PC3JyzNoKyjBT1IeaqrqczLI2HIYKujWEpESBo2DKLsQXsbtoolS6jMyyN90iRs8fGWZqnetp3iF1+kescOS3Moa2hBD1LOlSsx1dUkXHCB1VEs13ry/XSYOjVo76VqamqIy+1Hi2uvsToKqTf/ComLo3DmTKujKAtoQQ9SVRs2gMNBQgTOPz+dYL0aMumii+j8yitBMX4d1bIlLa8bx5H5H1K9c6fVcVSAaUEPUmm33ELXr/6LLSHB6ihB4eCTT7Jj9OVBNexi3G5K5871+eJbzZV2881IdDRF/5pldRQVYFrQg5gvl1wNddFZWdQWFFC9yfv7cfrbkfnz2f/Ag5R/9ZXVUU4QlZ5Oxt13k3TJMKujqADTgh6EyhYupOB3v4u4NdDPJPGii8Bmo2zhQqujAHXj5oeenUpsz551J22DTOqN40m65BKrY6gA04IehMoWfU7FsjzsyclWRwkaUampxPftS9lnwVHQS+fOpbaggFZ33en3Bbiayl1WxqFnp1JTUGB1FBUgXn0SRWSEiGwRkW0icspiESLSXUS+FZFqEbnX9zEjh/F4KP/qKxIHD0bsdqvjBJWk4ZdQvXWr5VeNepxOCqfPIC63HwlDhlia5Uw8lZUUzZpF0ax/Wx1FBchZF+cSETswDRgOFAB5IjLPGLPxuN2KgTuAK/wRMpJUrVuHu6iIxIuGWh0l6CRdeinY7JbfpcdVVERU27Zk3Hln0E6lhLqbYLS45mpK3nqb9EkTcbRrZ3Uk5Wfe9NAHANuMMTuMMTXA68AJq/YbYw4aY/KAWj9kjChlX34JNhuJP/mJ1VGCjqNtW1LH32D5XXqiMzPp/MbrxOfmWprDG2m33AJA0XPPW5xEBYI3Bb09kH/c44L65xpNRCaIyHIRWX7o0KGmHCLsOTIySLnqSsuLVrDyVFZSOncuNbt2WdJ+2Zdf4iouDuqe+fEc7drR4oorKH3rLWp/+MHqOMrPvCnoDX1ymzQZ2BgzyxiTa4zJbRXB65OcScvrrqPdo49aHSNoeZxO9j/0fymd+07A2649eJC9d97FwSeeDHjbzZE2cQIJQ4ZggvTCLOU73hT0AqDDcY8zAV102Q9qDxzAU1NjdYygFpWWRsKQwRyePx/j8QS07cIZMzAuF+m3Tgpou80VnZlJhxnTie7Uyeooys+8Keh5QFcRyRKRaGAcMM+/sSLT/oceYte1v7A6RtBLuXwMrv37A3oD6Zrduyl9621aXnsN0R07BqxdX6rJz+fIx59YHUP50VkLujHGBdwOfAJsAt40xmwQkUkiMglARNqISAFwN/CgiBSIiE6ibgT3kSNUfPsdCeefb3WUoJc07GJsycmUvvlGwNo89I9/Ig4H6bfeGrA2fa1w6lT2TZmi9x4NY17NQzfGLDDGdDPGnGOM+Uv9czONMTPrvz5gjMk0xiQbY1rUf33En8HDTdmiz6G2luSRI6yOEvRscXG0uPJKavcfwLjdfm/PuFwY4yHt5l+F9Nr0aRMnYqqqKH5xttVRlJ/oTaKDxJGPP8LRrh2x2dlWR2m2zpM/bPD5XY+N8lkbGffcjQRodUOJiiLz6acDPmbvazFdupA8ciQlc+aQ+uubda2gMBSc1yxHGPfhw1Qs+ZakESNCZjqc1Y4Wc3dpKcbl8ls75V8vpnrbtro2g/QS/8ZIv3USnspKiv/zH6ujKD8I/U9oGLAlJ9P5lZdpef11VkcJKVWbN7P1p0M58tHHfjm+u6yMfZMnc+DhP/nl+FaI6dqV5NGjMdU6myoc6ZBLEBAR4nr3tjpGyInp1o3ojh0omvUvkkdd5vMe9MGnnsJdUkLGv/7l0+Nard0Tf9O/BMOU9tAt5ioqYv9D/5eaPXusjhJyxGYjbcJEqrduo2zRIp8euzIvj9LXXif1xhuJ63WuT49ttaPFvHLVKtxlZRanUb6kPXSLHZk/n9I33yR1/A0+Pe7pTkyGm+SRIzg09VkK//ksSRddhEQ1/yPtqapi/4N/xJGZSas7fueDlMGnZvdudl93PWkTJ5Jx151Wx1E+oj10i5W++x6x2dnEdO1qdZSQJFFRZNx9D9XbtlG5YqWPDiokDb+Eto/8CVt8vG+OGWSiO3UiedQoil96idqDB62Oo3xEC7qFqjZtonrzZlKuvMLqKCEt6dLhdPlwPgkDBzT7WMYYbDExZNx7LwkXXOCDdMGr1e/vwLhcFE6bbnUU5SNa0C1U+u67iMNBymWXWR0lpIkIMV26AFC9dWuTj1OzZw87r7gS57r1vooW1KI7dqTltddS+vbbVO/caXUc5QNa0C1kT0wi5YordKlcHyn/+mt2XD6GIx83fhqjq6SE/AkTcR04gD0lclatSL/tVqIyMqjRgh4W9KSohcL1hJtVEgYNIjYnh/0PPEh0Vhdif9zNq9e5y8oomHQrtXv30nH2iyG7+FZTRKWn86PPPvXJyWRlPe2hW8AYQ+WqVRjTpGXl1WmIw0HmP/+BLSGB/AkTvBpGcB85wp5f3Yxz40baPfV34vv1C0DS4CJRURiPhyOffhqQtXGU/2hBt4Bz1Sp2X3c9R+ZHxtTCQHK0aUOHf8/C1Nay+8Yb8VRUnHF/W3w89pYtyXz2nyQPHx6glMGn/Kuv2HvH7yl9J/A3DlG+owXdAsWzX8KWnEzSsIutjhKWYn/8YzrNeYX0iZOwJSQAULFsGbU/HKR23z6OfPopBb+7A9ehQ0hUFB3+PYukoUOtDW2xxJ/+lLh+/Tj01NO4Dx+2Oo5qIh04C7DqHTso++wz0iZOCNs5zsEgJiuLmKwsACpXrmTPjTedsN2elkbVpk0ktmqll8FTN1OozYMPsPPnV3Po2am0efABqyOpJtCCHmBFs/6NxMSQeuONVkeJGLE9etDxhefrl1cQYn50DnF9+uiJwJPE9uhBy3G/oOTVV0kZczlxOTlWR1KNpJ/oAPJUV1OxbCktrr2GqNRUq+NEDFtcHAkXXBD2Fwr5Qqu776Zqy/cYvbdtSNKCHkC2mBjO+fhjTHW11VGUapA9MZHOc16xOoZqIj0pGiCuwkI8NTXYoqOxJyVZHUepM/JUV3Pw6Weo2rjR6iiqEbSgB8j+Bx5k97jrdO65CgnG6eTwu++y9+57zjr1UwUPHXIJgPJvvqH8v/8l4957fD6jIlKWyVWBZW/RgnZPPsGeX93MgUf+TLvHH7M6kvKC9tD9zON0cuDhPxHdqRMtx4+3Oo5SXksYMID0227j8PvvU/rue1bHUV7QHrqfFU6fQW1+Ph1nz8YWE2N1HEud7q+JXY+NCnAS5a30WydRuWwZBx97jKThw7EnJlgdSZ2BFnQ/Mh4PzjVrSLnyShIGDbQ6jlKNJnY77Z95mtr9+7WYhwAt6H4kNhsdX3hepymqkBaVmnrsuonD8z8k8cKfYE+OnCWGQ4mOofuBMYai51/AVVyMREUdW09EqVBWU1DAvilTyJ84SWe+BCkt6H5QOG06B594giPz51sdRSmfic7MpP2TT+Jcs4b8236rRT0IaUH3sdL33qNw6lRSrrhCZ7WosJP8s0tp99j/ozIvj903/xpXSYnVkdRxdAzdhw5/8AH7//AA8QMH0vaRP+kqfuq0QnnGT8qYMdgSEth79z1ULFlCyqjgzxwptKD7iKmpoXDGTOL796fDjOlIdLTVkZTym6Rhwzjnk49xtGkDQO3+/TjatrU4ldKC3kzG5cK43dhiYuj44ovYk5OwxcX5vB29IjQyhFLP/Wgxr9qyhV3XXEvK2LFk/M99ulaRhXQMvRmqd+xg1/W/ZP+Df8QYg6N1hl+KuVLBLDori5Y33EDp3LnsuHwMRz76SNcssoj20JvAXV5O0XPPUfzibGyxsaTd/CsdL1d+Fcw9d1t0NK3/5z6Sf3Yp+//4EHvvupu4OXPo9NJLiN1udbyIogW9kSqWLmPv3XfjLioiedQoMu7/HxwZGVbHCmnBXKyaK5KGyuJ69ybr3Xc4/P48XIWFiN2OMYayTz4lYcgQvdI0ALSgn4WnupqKxYuxt2hBfL9+RHfqSGyPHrS643fNukVXOBcxf9PvXfASu50WV1157HHVunXsvfNOJDqaxJ9eSNKlPyPh/EFEpadbmDJ8eVXQRWQE8A/ADjxnjHnspO1Sv/0yoBL4lTFmpY+zBkzFkiVUrlhJ5fLlOFevxlRXk3zZZcT364ejTRs6PvdvqyNGhGDs3QbbL5Ngy3Oy2OxsOr36Kkc++oiyjz+m7LOFAHT8z0skDBhATX4+roMHienWTU+m+sBZC7qI2IFpwHCgAMgTkXnGmONvZTIS6Fr/30BgRv3//coYA243xu0GtxtbfDxQN8btqawEjwfjcmOqqzC1tcR27w5AxbffUrNrF67CIlyHDuE6eBB7y5a0+39/BeDgk3+navNmYrt3p+W4X5DwkwtJGDjA328HCM4iFuoC8T3Vf7eGiQjxfc8jvu95tJ58P1UbNlCxdClx554LwOF58yh8dioAtuRkHG3b4mjXjnZP/A17YiKVK1ZQs2sXtuRk7ElJSEwMtpgYYnv2BMBdVgZuN0Q5EJuAzQY2G7b6acPGmIg6v+VND30AsM0YswNARF4HxgLHF/SxwH9M3ant70SkhYi0Ncbs93li4MCfH6Xk1VfhuDPpEhND9zWr67Y/8ghH5n1wwmvsqal0W/INAMWvzKF80aK659PSiMrIOGG9lfZP/R17WlqDPYbG/uAGS08pEmhR/V+B+Jw29q8DsduJy8k5Yaiy5S9+Qey551KzbRu1+/ZTu38/rh9+ONY5O/zBB5S+/saJx4mOpvvaNQD88OhfOPz++ydst7dsSbdvlwBQcPvv6n7W6wu9iODo2JFzPqxblmPPxIlULl12wutjunUj6826NnffMB7nhg0nbI/r05tOL74IwM6rfk71zp0nbE8YfAEdptb9kto+YiS1Bw+esD15+CW0e/zxBr9HzSVnm14kIlcDI4wxt9Q/Hg8MNMbcftw+84HHjDGL6x8vAu43xiw/6VgTgAn1D38MbPHVGwmgdKDQ6hABpu85/EXa+4XQfc+djDGtGtrgTQ+9ob9XTv4t4M0+GGNmAbO8aDNoichyY0yu1TkCSd9z+Iu09wvh+Z69ubCoAOhw3ONMYF8T9lFKKeVH3hT0PKCriGSJSDQwDph30j7zgBulziDgsL/Gz5VSSjXsrEMuxhiXiNwOfELdtMUXjDEbRGRS/faZwALqpixuo27a4s3+i2y5kB4yaiJ9z+Ev0t4vhOF7PutJUaWUUqFBF+dSSqkwoQVdKaXChBb0ZhCRe0XEiEhYL0whIk+IyGYRWSsi74pIC6sz+YuIjBCRLSKyTUQmW53H30Skg4h8ISKbRGSDiPze6kyBIiJ2EVlVfx1NWNCC3kQi0oG65RD2WJ0lAD4DehljcoDvgSkW5/GL45a5GAn0BK4TkZ7WpvI7F3CPMaYHMAj4bQS856N+D2yyOoQvaUFvuqeB/6GBC6jCjTHmU2OMq/7hd9RdZxCOji1zYYypAY4ucxG2jDH7jy6kZ4wpo67Atbc2lf+JSCYwCnjO6iy+pAW9CURkDLDXGLPG6iwW+DXwkdUh/KQ9kH/c4wIioLgdJSKdgfOApRZHCYRnqOuQeSzO4VO6HvppiMhCoE0Dmx4A/gBcGthE/nWm92uMeb9+nweo+xN9TiCzBZBXS1iEIxFJBOYCdxpjjlidx59EZDRw0BizQkSGWhzHp7Sgn4Yx5pKGnheRbCALWFO/LGcmsFJEBhhjDgQwok+d7v0eJSI3AaOBYSZ8L16IyCUsRMRBXTGfY4x5x+o8ATAYGCMilwGxQLKIvGKMucHiXM2mFxY1k4jsAnKNMaG4aptX6m9w8hTwU2PMIavz+IuIRFF30ncYsJe6ZS+uN8ZsOOMLQ1j9zWleAoqNMXdaHCfg6nvo9xpjRlscxSd0DF15YyqQBHwmIqtFZKbVgfyh/sTv0WUuNgFvhnMxrzcYGA9cXP9vu7q+56pCkPbQlVIqTGgPXSmlwoQWdKWUChNa0JVSKkxoQVdKqTChBV0ppcKEFnSllAoTWtCVUipM/H9eI9kt0jI/rAAAAABJRU5ErkJggg==",
      "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",
	" self.t = 1.0\n",
	" self.rho = 0.999\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) ** (1.0/self.t)\n",
	" p1 = self.inv_dist_fun(x1) ** (1.0/self.t)\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",
	" self.t *= self.rho\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/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtyklEQVR4nO3deXxU1f3/8ddnJpN9gSSELSzBQgFJQAiLQi2KWBAEtWrRitbWH6C11vUrVOvPWtufVqu2spW6YBV3XBBxA7WKKIR9R3YSFiEbZJksM3N+fyTwZQkwSWbmzvJ5Ph4+zMy9c897wuSTk3PPPVeMMSillAp9NqsDKKWU8g0t6EopFSa0oCulVJjQgq6UUmFCC7pSSoUJLehKKRUmvCroIjJCRLaIyDYRmdzA9qEiclhEVtf/95DvoyqllDqTqLPtICJ2YBowHCgA8kRknjFm40m7fm2MGe2HjEoppbzgTQ99ALDNGLPDGFMDvA6M9W8spZRSjXXWHjrQHsg/7nEBMLCB/c4XkTXAPuBeY8yGk3cQkQnABICEhIR+3bt3b3xiFXDr9h5u8Pns9ikBTqKUWrFiRaExplVD27wp6NLAcyevF7AS6GSMKReRy4D3gK6nvMiYWcAsgNzcXLN8+XIvmldW6zz5wwafX/7YqAAnUUqJyO7TbfNmyKUA6HDc40zqeuHHGGOOGGPK679eADhEJL0JWZVSSjWRNwU9D+gqIlkiEg2MA+Ydv4OItBERqf96QP1xi3wdViml1OmddcjFGOMSkduBTwA78IIxZoOITKrfPhO4GrhVRFyAExhndBlHpZQKKG/G0I8Ooyw46bmZx309FZjq22hKqXBSW1tLQUEBVVVVVkcJCbGxsWRmZuJwOLx+jVcFXSmlmqugoICkpCQ6d+5M/QitOg1jDEVFRRQUFJCVleX16/TSf6VUQFRVVZGWlqbF3AsiQlpaWqP/mtGCrpQKGC3m3mvK90oLulJKhQkt6EqpiFBaWsr06dP93s57773Hxo0nL3UVGFrQlVIRobEF3RiDx+NpdDta0JVSys8mT57M9u3b6dOnD3fddRfDhg2jb9++ZGdn8/777wOwa9cuevTowW233Ubfvn3Jz8/nz3/+M927d2f48OFcd911PPnkkwBs376dESNG0K9fP37yk5+wefNmlixZwrx587jvvvvo06cP27dvD+h71GmLSilL7B5/4ynPJY0cQer11+NxOsmfMPGU7SlXXkmLq67EVVLC3jt+f8K2Ti//54ztPfbYY6xfv57Vq1fjcrmorKwkOTmZwsJCBg0axJgxYwDYsmULL774ItOnT2f58uXMnTuXVatW4XK56Nu3L/369QNgwoQJzJw5k65du7J06VJuu+02Pv/8c8aMGcPo0aO5+uqrm/qtaTIt6EqpiGOM4Q9/+ANfffUVNpuNvXv38sMPPwDQqVMnBg0aBMDixYsZO3YscXFxAFx++eUAlJeXs2TJEq655ppjx6yurg7wuziVFnSllCXO1KO2xcWdcXtUy5Zn7ZGfyZw5czh06BArVqzA4XDQuXPnY3O+ExISju13uhVMPB4PLVq0YPXq1U3O4A86hq6UighJSUmUlZUBcPjwYTIyMnA4HHzxxRfs3t3wirRDhgzhgw8+oKqqivLycj78sG4p6eTkZLKysnjrrbeAusK/Zs2aU9oJNC3oSqmIkJaWxuDBg+nVqxerV69m+fLl5ObmMmfOHE53s53+/fszZswYevfuzVVXXUVubi4pKXU3dpkzZw7PP/88vXv35txzzz12YnXcuHE88cQTnHfeeQE/KSpWLYqoN7gIHae7wcUuvcGFaoRNmzbRo0cPq2M0Wnl5OYmJiVRWVnLhhRcya9Ys+vbtG5C2G/qeicgKY0xuQ/vrGLpSSp3BhAkT2LhxI1VVVdx0000BK+ZNoQVdKaXO4NVXX7U6gtd0DF0ppcKEFnSllAoTWtCVUipMaEFXSqkwoSdFlVKWON102KYKpmm0Q4cO5cknnyQ3t8HZhX6jPXSllAoTWtCVUhGhoqKCUaNG0bt3b3r16sUbb7zBI488Qv/+/enVqxcTJkw4tnbL0KFDueuuu7jwwgvp0aMHeXl5XHXVVXTt2pUHH3wQqFtqt3v37tx0003k5ORw9dVXU1lZeUq7n376Keeffz59+/blmmuuoby8HKhbzrdnz57k5ORw7733+uQ9akFXSkWEjz/+mHbt2rFmzRrWr1/PiBEjuP3228nLy2P9+vU4nU7mz59/bP/o6Gi++uorJk2axNixY5k2bRrr169n9uzZFBUVAXVL7U6YMIG1a9eSnJx8yg00CgsLefTRR1m4cCErV64kNzeXp556iuLiYt599102bNjA2rVrj/2SaC4t6EqpiJCdnc3ChQu5//77+frrr0lJSeGLL75g4MCBZGdn8/nnn7Nhw4Zj+x9dHz07O5tzzz2Xtm3bEhMTQ5cuXcjPzwegQ4cODB48GIAbbriBxYsXn9Dmd999x8aNGxk8eDB9+vThpZdeYvfu3SQnJxMbG8stt9zCO++8Q3x8vE/eo54UVcf4+iSVUsGkW7durFixggULFjBlyhQuvfRSpk2bxvLly+nQoQMPP/zwsSV0AWJiYgCw2WzHvj762OVyASAiJ7Rx8mNjDMOHD+e11147Jc+yZctYtGgRr7/+OlOnTuXzzz9v9nvUHrpSKiLs27eP+Ph4brjhBu69915WrlwJQHp6OuXl5bz99tuNPuaePXv49ttvAXjttdcYMmTICdsHDRrEN998w7Zt2wCorKzk+++/p7y8nMOHD3PZZZfxzDPP+Gxdde2hK6UsEehphuvWreO+++7DZrPhcDiYMWMG7733HtnZ2XTu3Jn+/fs3+pg9evTgpZdeYuLEiXTt2pVbb731hO2tWrVi9uzZXHfddcfuaPToo4+SlJTE2LFjqaqqwhjD008/7ZP3qMvnqmMaO+QSTPN+VfAL1eVzT2fXrl2MHj2a9evX+62Nxi6fq0MuSikVJrSgK6VUE3Tu3NmvvfOm0IKulAoYq4Z4Q1FTvlda0JVSAREbG0tRUZEWdS8YYygqKiI2NrZRr9NZLkqpgMjMzKSgoIBDhw5ZHSUkxMbGkpmZ2ajXaEFXSgWEw+EgKyvL6hhhTYdclFIqTHhV0EVkhIhsEZFtIjL5DPv1FxG3iFztu4hKKaW8cdaCLiJ2YBowEugJXCciPU+z3+PAJ74OqZRS6uy86aEPALYZY3YYY2qA14GxDez3O2AucNCH+ZRSSnnJm4LeHsg/7nFB/XPHiEh74Epg5pkOJCITRGS5iCzXM91KKeVb3hR0aeC5kyeSPgPcb4xxn+lAxphZxphcY0xuq1atvIyorBRXW8Vv18yle/Fuq6Mopc7Cm2mLBUCH4x5nAvtO2icXeL1+LeB04DIRcRlj3vNFSBVYXUr30qaymCXtsmlbUcTF+SsZvfNb3ux6EbN7jsSITo5SKhh5U9DzgK4ikgXsBcYB1x+/gzHm2ORSEZkNzNdiHpralhfy1yWzKIuOZ2mbnuxo0Z5fjniI/7P+A67d+gUesfFSz5FWx1RKNeCsBd0Y4xKR26mbvWIHXjDGbBCRSfXbzzhurkKH3eNmSt7LgOGhQb/BbbMDUBUVw7O9f44Yw7jvF7EptRPL2pwy0UkpZTGvrhQ1xiwAFpz0XIOF3Bjzq+bHUlb4+bYv6Xp4L3/pP579ieknbhRheu8rAchPzLAgnVLqbPTSfwWAq6iIcVsWsaRtLxa3793wPrYo/nneNQFOppTylhZ0BYC9RQum9b6KzamdzrpvRkUx4zd/gqv4fKJSUwOQTinlDZ2uoAAQu51FHXPZm3j26aQxnlouzl9JySuvBCCZUspbWtAVJW+9ReGsf4OX61TnJ7Umr3V3Sl5/A0/9jW+VUtbTgh7hjNtN4YwZVHzzDUhD15A17L0fXYi7uJgj8xt3Y2mllP9oQY9w5V9+iWvfflpef/3Zdz7O6vQfEX3OOZS+/bafkimlGksLeoQrnfsOUa1akTTs4sa9UITU8eOJy8nGuM+44oNSKkB0lksEcx8+TPnXX5P6y18iUY3/KLQc9ws/pFJKNZX20COY+/BhEgcPJnn06CYfw3g8VOblYTweHyZTSjWFFvQIFt2xIx1mziCu17lNPsaRjz5i9/gbqVq3zofJlFJNoQU9QrnLy6k9cKDZx0kcMgSioihbuNAHqZRSzaEFPUKVffIJ24ZeRPWOHc06jj0lhYQBAyj7bCHGy3nsSin/0IIeoco+/4Kotm2Jzso6+85nkTT8Emp27aJm+3YfJFNKNZUW9Ajkqa6mYskSki4aijTiYqLTSbx4GABlCxc1+1hKqabTaYsRqPK77zBOJ4kXNXLu+Wk4WmfQ+a03if3xj31yPKVU02hBj0BlX3yBLT6e+IEDfHbMuOxsnx1LKdU0OuQSgdJvu432z/4TW3S0z47pLq/g4JNPUv7NNz47plKqcbSHHoEcGRk4Mnx71yFbbAwlb76Fq6SExMGDfXpspZR3tIceYcq//priOXN8vv6KREWRcMEFVHy9WKcvKmURLegRpuSNNyh+4UWw+f6fPnHIYFwHD1K9davPj62UOjst6BHEuN1ULl1G/PmDfDJd8WQJQ4YAULFYx9GVsoIW9AhStWEDnrIyEs4/3y/Hd7RpQ1zv3piaGr8cXyl1ZnpSNIJUfPsdAAmDBvmtjc5vvO63Yyulzkx76BGkdv8+Ynr0ICotze9t6YlRpQJPC3oEafvww2T5uQftqaxk+2WjKJ79kl/bUUqdSgt6hBEfXkzUEFt8PHg8VC5b5td2lFKn0oIeIYr/8zJ7bvk/mNpav7cVP2AAlcuX671GlQowLegRonzx19Qe2I84HH5vK37AADxlZVRt3uz3tpRS/0sLegQwbjfOVauJ79svIO3F9+8PQOVSHXZRKpC0oEeA6m3b8JSVEZ8bmILuaJ1B6m9+TcyPuwWkPaVUHZ2HHgEqly8HIC5APXSA1vfdF7C2lFJ1tIceAaLS0kj62c9wtG8XsDaNMVTv3In78OGAtalUpNOCHgGSR4wg8x/P+GX9ltOp3rqVHSMvo+yLLwLWplKRTgt6mPM4nXgqKgLebsyPfoQtMRHnqtUBb1upSOVVQReRESKyRUS2icjkBraPFZG1IrJaRJaLyBDfR1VNUfbpp2wZMJCaXbsC2q7YbMT16YNz1aqAtqtUJDtrQRcROzANGAn0BK4TkZ4n7bYI6G2M6QP8GnjOxzlVE1WuXo0tNhZHhw4BbzvuvD5Ub92Ku6ws4G0rFYm86aEPALYZY3YYY2qA14Gxx+9gjCk3/7saUwKgKzMFiaq164jNzkbs9oC3Hd+3LxiDc83agLetVCTypqC3B/KPe1xQ/9wJRORKEdkMfEhdL/0UIjKhfkhm+aFDh5qSVzWCp7qaqi1biMvOtqT9uN69af/sP4nL7mVJ+0pFGm8KekNTI07pgRtj3jXGdAeuAP7c0IGMMbOMMbnGmNxWrVo1KqhqvOpNm8DlIjbHmoJui48nefhw7CkplrSvVKTxpqAXAMcPwGYC+063szHmK+AcEUlvZjbVTFFt29J6yuS6oQ+L1OzeTdGLs3WhLqUCwJuCngd0FZEsEYkGxgHzjt9BRH4k9ZOcRaQvEA0U+TqsahxH69ak3nRTQG5ocTrOtWs5+PjjeuNopQLgrAXdGOMCbgc+ATYBbxpjNojIJBGZVL/bz4H1IrKauhkxvzB6yxrLlf/3v7gsPlcRd955ADp9UakA8GoeujFmgTGmmzHmHGPMX+qfm2mMmVn/9ePGmHONMX2MMecbYxb7M7Q6O3dpKfkTJ1H6zruW5nC0b489NRXnuvWW5lAqEuiVomHKuX4DAHEWnRA9SkSIy87GuXaNpTmUigRa0MNU1bq6ud+xvayfMhibk03t7j14nE6roygV1rSghynn2nVEd+mCPSnJ6iik3ngj3fKWYYuLszqKUmFNC3oYMsbgXLfOsguKTmZPSsIWG2t1DKXCnt7gIkx1fnUOeDxWxzimaPZsPEfKaHXH76yOolTY0h56GBIRojt2JLpzZ6ujHFO1fgOlb79tdQylwpoW9DB05KOPgq54xuXk4Dp4kNoffrA6ilJhSwt6GCp59TVK3nzL6hgnODp90rlWV15Uyl+0oIcZ43bj3LCBuJwcq6OcIKZHD3A4qNKCrpTfaEEPM9XbtmMqKy2/oOhktpgYEgYMsDqGUmFNZ7mEmWMXFAXJlMXjdXxeb2SllD9pDz3M1BQUYEtJIbpTJ6ujKKUCTAt6mMm48066fvkFYgu+f1pXYSE7Lr+c0nffszqKUmEp+H7qVbMF6yX29tRUag/8gHP1aqujKBWWtKCHEee69eyZOJHqnTutjtIgsdmIy+6Fc53OdFHKH7SghxHnyhVU/PcrbAkJVkc5rdjsHKq3fK8rLyrlB1rQw4hz7Tqi2rbFkZFhdZTTiuudA243VRs3Wh1FqbCjBT2MONeuDZoVFk8nLieH5DGXY4uPtzqKUmFHC3qYcJWUUJufH3QXFJ0sKj2d9n/7G7E9elgdRamwowU9TLhLSonPzT12U+ZgZoyhdt8+q2MoFXa0oIeJmC5ZdHrlZeL79bM6ylmVvPwK2y4ehqu42OooSoUVLehhwrjdVkfwWmzPuuEWXXlRKd/Sgh4GjDFsu+hiDk2dZnUUr8T27Ak2m668qJSPaUEPA7V79+I6eJCo9DSro3jFFh9PTNeuONeuszqKUmFFC3oYONrTDcYVFk8nLicH57p1GGOsjqJU2NDlc8OAc+06JCaG2G7drI7itRZX/5yEC84Htxui9GOolC/oT1IYcK5bR2zPnojDYXUUr8X17k1c795Wx1AqrGhBDwPJI0YE9fotp+NcvwFPZYXeyUgpH9GCHgZSx99gdYQmOfj443iqq8l68w2roygVFvSkaIir3bcPV2Gh1TGaJDYnm+pNm/DU1FgdRamwoAU9xB2aPp0do0aH5GyRuJzemNpaqjdvtjqKUmFBC3qIq1q7jtjsbETE6iiNdnQhMecavcBIKV/QMfQQ5qmooHrbNpIuucSS9jtP/rDB53c9Nsqr10e1aUNUq1a6BIBSPqIFPYQ5128Aj6fuphEhSETo8PxzONq1tzqKUmFBC3oIc65dA0BsTmgWdCCkLoZSKth5VdBFZATwD8AOPGeMeeyk7b8E7q9/WA7caoxZ48ug6lQpl11GdGYmUS1bNup1pxsqsYKrpISSl18m8aKLicvuZXUcpULaWQu6iNiBacBwoADIE5F5xpjjbwq5E/ipMaZEREYCs4CB/gis/pejfXsc7UN7uEKioiicMROiorSgK9VM3sxyGQBsM8bsMMbUAK8DY4/fwRizxBhTUv/wOyDTtzHVyVwlJZS+/TauoiKrozSLPSmJ6C5dqNKVF5VqNm8Kensg/7jHBfXPnc5vgI8a2iAiE0RkuYgsP3TokPcp1Skqly9n/4N/pDY//+w7B7m4nByca9eG5Fx6pYKJNwW9oQnODf7kichF1BX0+xvaboyZZYzJNcbktmrVyvuU6hRVa9eCw0FMGNxsOS4nG3dxMbV79T6jSjWHNwW9AOhw3ONM4JSfPBHJAZ4DxhpjQnscIAQ416wltnt3bDExVkdptticHGzx8dQWhP5fG0pZyZtZLnlAVxHJAvYC44Drj99BRDoC7wDjjTHf+zylOoFxu6lav56UK66wOopPxPboQbe8ZYjdbnUUpULaWQu6McYlIrcDn1A3bfEFY8wGEZlUv30m8BCQBkyvvwTdZYzJ9V/syFazezeeysqQvaDoZGLTFSiU8gWv5qEbYxYAC056buZxX98C3OLbaOp0Yrp0oeuSb5Do0B9uOerIZ59RNOvfdJrzCrboaKvjKBWStGsUoqJSU7Enht5NLc6kat06qtZvsDqGUiFLC3oIOvDXv3LkowZnhoas+PPOA8C5aqXFSZQKXVrQQ4ynspKSV+ZQvXWr1VF8Kio9HUenjlSuXGV1FKVClhb0EFO1oW6FxVBekOt04s/ri3PlSr3ASKkm0tUWQ4xzTd2aZ3G9e1ucxPcSf3ohpqYaT0Vl2J0fUCoQtKCHmMqVq4ju3LnRKyyGguSRI0keOdLqGEqFLC3oocYY4geF90KW7rIy7ElJVscImJqCvVQs+YbaPXtIHDbs2AlipRpLC3qI6TBjeliPMe+7fzLOtWs556MFZ985xNXs3s0Pf3uC8s8/B2MQh4P4gYMAqFy5kqLnXyDj7ruIOecci5OqUKEnRUNQKN4Q2lvRXbpQs3MnrpKSs+8cwtzlFey85loqly4l/dZJdFmwgB+vXUPiT4YAULt3H5XLlrHzyqsofmVOWP8SV76jBT2EHHjkz+RPutXqGH4V3/fofPTwnr5oT0yg7aN/psuHH9LqjjuI6ZJ1wi/qlMtHc85HC0i44AJ+ePRRDj72GMbjsTCxCgVa0ENIxdKlEOY9tdjsbHA4cK4MvwuMjDEc+OtfOTy/7haAyZdeiqN1xmn3j0pPJ3P6NFJvupHil/7D4XnzAhVVhSgt6CHCXVpKzfbtxPXta3UUv7LFxhLbs0dYXmBUOHUaJf95merNm7x+jdhstJ4yhcyZM0gZM8aP6VQ40JOiIaJy9WoA4s7rY2mOQEj7zW/AE15/iZQtWkThtGmkXHEFre65p9GvTxo6FKibEVO7dy8JAwf4OKEKB1rQQ4Rz5ar6GylnWx3F75IvvdTqCD5Vs2sX++6fTOy559LmTw8366T2/j8+SNXGTWTNnUt0ZmjfIFz5ng65hIiYbt1I/eUvscXFWR0lIKq2fB82wy5liz5HoqLI/Oc/mn2HqbZ/+hN4POy75x6M2+2jhCpcaEEPESmjR9F6ymSrYwTM/of+yMG//c3qGD6R9ptf0+WjBTjaN79HHd2xI20e+iPONWsofvllH6RT4UQLeghwlZSE/bzskyUMGIhz/Xo8FRVWR2my6h07cK5bD+DTpRqSR48mcehQDj3zD2p27/bZcVXo04IeAkrfeputFwzGXVpqdZSAiR84EFyukB12MR4P+//wAPm33YqnutqnxxYR2vzpYVr+4lrsKSk+PbYKbVrQQ0BlXh7R53TB3qKF1VECJr7veeBwULlsqdVRmuTw+/Nwrl5Nxl13N3vcvCGO1q1pPWVKRH0m1NlpQQ9yxuXCuWIFCQMia5qaLT6euOxsKpflWR2l0dxlZRz8+9+J692blCvG+rUt57p15N9+O56qKr+2o0KDFvQgV7VxI57KSuL797c6SsC1/cujdHj+OatjNFrhtOm4i4po/eCDiM2/P2Ke8nLKFy6i+MUX/dqOCg1a0INc5bJlABFZ0GOysrAnJlodo9HsLVvS8oYbiMvu5fe2Es4/n8RLhlH03PO4iov93p4KbnphUZBLGjGSqNatiUpPtzqKJYpmz0ZsdlJvHG91FK+lT5wQ0PYy7rqLHZ+PoehfsyJqaqs6lfbQg1x0ZntSLr/c6hiWqfz2O4rnvGJ1DK9Ub91K2cKFAV/qNuacc0i58gpKXnuN2n37Atq2Ci5a0INYzZ49lL77Hu7ycqujWCZhyBBqd+8JifnWB//+FPv+8AAeC/69Wv32t2Tcczf2CP1LTtXRgh7Eyj5byP4pUzBOp9VRLHP0hg/lixdbnOTMKleuovzLL0m75RZLbp/naNeO1JtuwhYdHfC2VfDQgh7EKr79luguXYhq1crqKJZxdOqEo0MHKr4O3oJujOHQ009jT08n9YZfWprl8Lx5/PC3JyzNoKyjBT1IeaqrqczLI2HIYKujWEpESBo2DKLsQXsbtoolS6jMyyN90iRs8fGWZqnetp3iF1+kescOS3Moa2hBD1LOlSsx1dUkXHCB1VEs13ry/XSYOjVo76VqamqIy+1Hi2uvsToKqTf/ComLo3DmTKujKAtoQQ9SVRs2gMNBQgTOPz+dYL0aMumii+j8yitBMX4d1bIlLa8bx5H5H1K9c6fVcVSAaUEPUmm33ELXr/6LLSHB6ihB4eCTT7Jj9OVBNexi3G5K5871+eJbzZV2881IdDRF/5pldRQVYFrQg5gvl1wNddFZWdQWFFC9yfv7cfrbkfnz2f/Ag5R/9ZXVUU4QlZ5Oxt13k3TJMKujqADTgh6EyhYupOB3v4u4NdDPJPGii8Bmo2zhQqujAHXj5oeenUpsz551J22DTOqN40m65BKrY6gA04IehMoWfU7FsjzsyclWRwkaUampxPftS9lnwVHQS+fOpbaggFZ33en3Bbiayl1WxqFnp1JTUGB1FBUgXn0SRWSEiGwRkW0icspiESLSXUS+FZFqEbnX9zEjh/F4KP/qKxIHD0bsdqvjBJWk4ZdQvXWr5VeNepxOCqfPIC63HwlDhlia5Uw8lZUUzZpF0ax/Wx1FBchZF+cSETswDRgOFAB5IjLPGLPxuN2KgTuAK/wRMpJUrVuHu6iIxIuGWh0l6CRdeinY7JbfpcdVVERU27Zk3Hln0E6lhLqbYLS45mpK3nqb9EkTcbRrZ3Uk5Wfe9NAHANuMMTuMMTXA68AJq/YbYw4aY/KAWj9kjChlX34JNhuJP/mJ1VGCjqNtW1LH32D5XXqiMzPp/MbrxOfmWprDG2m33AJA0XPPW5xEBYI3Bb09kH/c44L65xpNRCaIyHIRWX7o0KGmHCLsOTIySLnqSsuLVrDyVFZSOncuNbt2WdJ+2Zdf4iouDuqe+fEc7drR4oorKH3rLWp/+MHqOMrPvCnoDX1ymzQZ2BgzyxiTa4zJbRXB65OcScvrrqPdo49aHSNoeZxO9j/0fymd+07A2649eJC9d97FwSeeDHjbzZE2cQIJQ4ZggvTCLOU73hT0AqDDcY8zAV102Q9qDxzAU1NjdYygFpWWRsKQwRyePx/j8QS07cIZMzAuF+m3Tgpou80VnZlJhxnTie7Uyeooys+8Keh5QFcRyRKRaGAcMM+/sSLT/oceYte1v7A6RtBLuXwMrv37A3oD6Zrduyl9621aXnsN0R07BqxdX6rJz+fIx59YHUP50VkLujHGBdwOfAJsAt40xmwQkUkiMglARNqISAFwN/CgiBSIiE6ibgT3kSNUfPsdCeefb3WUoJc07GJsycmUvvlGwNo89I9/Ig4H6bfeGrA2fa1w6lT2TZmi9x4NY17NQzfGLDDGdDPGnGOM+Uv9czONMTPrvz5gjMk0xiQbY1rUf33En8HDTdmiz6G2luSRI6yOEvRscXG0uPJKavcfwLjdfm/PuFwY4yHt5l+F9Nr0aRMnYqqqKH5xttVRlJ/oTaKDxJGPP8LRrh2x2dlWR2m2zpM/bPD5XY+N8lkbGffcjQRodUOJiiLz6acDPmbvazFdupA8ciQlc+aQ+uubda2gMBSc1yxHGPfhw1Qs+ZakESNCZjqc1Y4Wc3dpKcbl8ls75V8vpnrbtro2g/QS/8ZIv3USnspKiv/zH6ujKD8I/U9oGLAlJ9P5lZdpef11VkcJKVWbN7P1p0M58tHHfjm+u6yMfZMnc+DhP/nl+FaI6dqV5NGjMdU6myoc6ZBLEBAR4nr3tjpGyInp1o3ojh0omvUvkkdd5vMe9MGnnsJdUkLGv/7l0+Nard0Tf9O/BMOU9tAt5ioqYv9D/5eaPXusjhJyxGYjbcJEqrduo2zRIp8euzIvj9LXXif1xhuJ63WuT49ttaPFvHLVKtxlZRanUb6kPXSLHZk/n9I33yR1/A0+Pe7pTkyGm+SRIzg09VkK//ksSRddhEQ1/yPtqapi/4N/xJGZSas7fueDlMGnZvdudl93PWkTJ5Jx151Wx1E+oj10i5W++x6x2dnEdO1qdZSQJFFRZNx9D9XbtlG5YqWPDiokDb+Eto/8CVt8vG+OGWSiO3UiedQoil96idqDB62Oo3xEC7qFqjZtonrzZlKuvMLqKCEt6dLhdPlwPgkDBzT7WMYYbDExZNx7LwkXXOCDdMGr1e/vwLhcFE6bbnUU5SNa0C1U+u67iMNBymWXWR0lpIkIMV26AFC9dWuTj1OzZw87r7gS57r1vooW1KI7dqTltddS+vbbVO/caXUc5QNa0C1kT0wi5YordKlcHyn/+mt2XD6GIx83fhqjq6SE/AkTcR04gD0lclatSL/tVqIyMqjRgh4W9KSohcL1hJtVEgYNIjYnh/0PPEh0Vhdif9zNq9e5y8oomHQrtXv30nH2iyG7+FZTRKWn86PPPvXJyWRlPe2hW8AYQ+WqVRjTpGXl1WmIw0HmP/+BLSGB/AkTvBpGcB85wp5f3Yxz40baPfV34vv1C0DS4CJRURiPhyOffhqQtXGU/2hBt4Bz1Sp2X3c9R+ZHxtTCQHK0aUOHf8/C1Nay+8Yb8VRUnHF/W3w89pYtyXz2nyQPHx6glMGn/Kuv2HvH7yl9J/A3DlG+owXdAsWzX8KWnEzSsIutjhKWYn/8YzrNeYX0iZOwJSQAULFsGbU/HKR23z6OfPopBb+7A9ehQ0hUFB3+PYukoUOtDW2xxJ/+lLh+/Tj01NO4Dx+2Oo5qIh04C7DqHTso++wz0iZOCNs5zsEgJiuLmKwsACpXrmTPjTedsN2elkbVpk0ktmqll8FTN1OozYMPsPPnV3Po2am0efABqyOpJtCCHmBFs/6NxMSQeuONVkeJGLE9etDxhefrl1cQYn50DnF9+uiJwJPE9uhBy3G/oOTVV0kZczlxOTlWR1KNpJ/oAPJUV1OxbCktrr2GqNRUq+NEDFtcHAkXXBD2Fwr5Qqu776Zqy/cYvbdtSNKCHkC2mBjO+fhjTHW11VGUapA9MZHOc16xOoZqIj0pGiCuwkI8NTXYoqOxJyVZHUepM/JUV3Pw6Weo2rjR6iiqEbSgB8j+Bx5k97jrdO65CgnG6eTwu++y9+57zjr1UwUPHXIJgPJvvqH8v/8l4957fD6jIlKWyVWBZW/RgnZPPsGeX93MgUf+TLvHH7M6kvKC9tD9zON0cuDhPxHdqRMtx4+3Oo5SXksYMID0227j8PvvU/rue1bHUV7QHrqfFU6fQW1+Ph1nz8YWE2N1HEud7q+JXY+NCnAS5a30WydRuWwZBx97jKThw7EnJlgdSZ2BFnQ/Mh4PzjVrSLnyShIGDbQ6jlKNJnY77Z95mtr9+7WYhwAt6H4kNhsdX3hepymqkBaVmnrsuonD8z8k8cKfYE+OnCWGQ4mOofuBMYai51/AVVyMREUdW09EqVBWU1DAvilTyJ84SWe+BCkt6H5QOG06B594giPz51sdRSmfic7MpP2TT+Jcs4b8236rRT0IaUH3sdL33qNw6lRSrrhCZ7WosJP8s0tp99j/ozIvj903/xpXSYnVkdRxdAzdhw5/8AH7//AA8QMH0vaRP+kqfuq0QnnGT8qYMdgSEth79z1ULFlCyqjgzxwptKD7iKmpoXDGTOL796fDjOlIdLTVkZTym6Rhwzjnk49xtGkDQO3+/TjatrU4ldKC3kzG5cK43dhiYuj44ovYk5OwxcX5vB29IjQyhFLP/Wgxr9qyhV3XXEvK2LFk/M99ulaRhXQMvRmqd+xg1/W/ZP+Df8QYg6N1hl+KuVLBLDori5Y33EDp3LnsuHwMRz76SNcssoj20JvAXV5O0XPPUfzibGyxsaTd/CsdL1d+Fcw9d1t0NK3/5z6Sf3Yp+//4EHvvupu4OXPo9NJLiN1udbyIogW9kSqWLmPv3XfjLioiedQoMu7/HxwZGVbHCmnBXKyaK5KGyuJ69ybr3Xc4/P48XIWFiN2OMYayTz4lYcgQvdI0ALSgn4WnupqKxYuxt2hBfL9+RHfqSGyPHrS643fNukVXOBcxf9PvXfASu50WV1157HHVunXsvfNOJDqaxJ9eSNKlPyPh/EFEpadbmDJ8eVXQRWQE8A/ADjxnjHnspO1Sv/0yoBL4lTFmpY+zBkzFkiVUrlhJ5fLlOFevxlRXk3zZZcT364ejTRs6PvdvqyNGhGDs3QbbL5Ngy3Oy2OxsOr36Kkc++oiyjz+m7LOFAHT8z0skDBhATX4+roMHienWTU+m+sBZC7qI2IFpwHCgAMgTkXnGmONvZTIS6Fr/30BgRv3//coYA243xu0GtxtbfDxQN8btqawEjwfjcmOqqzC1tcR27w5AxbffUrNrF67CIlyHDuE6eBB7y5a0+39/BeDgk3+navNmYrt3p+W4X5DwkwtJGDjA328HCM4iFuoC8T3Vf7eGiQjxfc8jvu95tJ58P1UbNlCxdClx554LwOF58yh8dioAtuRkHG3b4mjXjnZP/A17YiKVK1ZQs2sXtuRk7ElJSEwMtpgYYnv2BMBdVgZuN0Q5EJuAzQY2G7b6acPGmIg6v+VND30AsM0YswNARF4HxgLHF/SxwH9M3ant70SkhYi0Ncbs93li4MCfH6Xk1VfhuDPpEhND9zWr67Y/8ghH5n1wwmvsqal0W/INAMWvzKF80aK659PSiMrIOGG9lfZP/R17WlqDPYbG/uAGS08pEmhR/V+B+Jw29q8DsduJy8k5Yaiy5S9+Qey551KzbRu1+/ZTu38/rh9+ONY5O/zBB5S+/saJx4mOpvvaNQD88OhfOPz++ydst7dsSbdvlwBQcPvv6n7W6wu9iODo2JFzPqxblmPPxIlULl12wutjunUj6826NnffMB7nhg0nbI/r05tOL74IwM6rfk71zp0nbE8YfAEdptb9kto+YiS1Bw+esD15+CW0e/zxBr9HzSVnm14kIlcDI4wxt9Q/Hg8MNMbcftw+84HHjDGL6x8vAu43xiw/6VgTgAn1D38MbPHVGwmgdKDQ6hABpu85/EXa+4XQfc+djDGtGtrgTQ+9ob9XTv4t4M0+GGNmAbO8aDNoichyY0yu1TkCSd9z+Iu09wvh+Z69ubCoAOhw3ONMYF8T9lFKKeVH3hT0PKCriGSJSDQwDph30j7zgBulziDgsL/Gz5VSSjXsrEMuxhiXiNwOfELdtMUXjDEbRGRS/faZwALqpixuo27a4s3+i2y5kB4yaiJ9z+Ev0t4vhOF7PutJUaWUUqFBF+dSSqkwoQVdKaXChBb0ZhCRe0XEiEhYL0whIk+IyGYRWSsi74pIC6sz+YuIjBCRLSKyTUQmW53H30Skg4h8ISKbRGSDiPze6kyBIiJ2EVlVfx1NWNCC3kQi0oG65RD2WJ0lAD4DehljcoDvgSkW5/GL45a5GAn0BK4TkZ7WpvI7F3CPMaYHMAj4bQS856N+D2yyOoQvaUFvuqeB/6GBC6jCjTHmU2OMq/7hd9RdZxCOji1zYYypAY4ucxG2jDH7jy6kZ4wpo67Atbc2lf+JSCYwCnjO6iy+pAW9CURkDLDXGLPG6iwW+DXwkdUh/KQ9kH/c4wIioLgdJSKdgfOApRZHCYRnqOuQeSzO4VO6HvppiMhCoE0Dmx4A/gBcGthE/nWm92uMeb9+nweo+xN9TiCzBZBXS1iEIxFJBOYCdxpjjlidx59EZDRw0BizQkSGWhzHp7Sgn4Yx5pKGnheRbCALWFO/LGcmsFJEBhhjDgQwok+d7v0eJSI3AaOBYSZ8L16IyCUsRMRBXTGfY4x5x+o8ATAYGCMilwGxQLKIvGKMucHiXM2mFxY1k4jsAnKNMaG4aptX6m9w8hTwU2PMIavz+IuIRFF30ncYsJe6ZS+uN8ZsOOMLQ1j9zWleAoqNMXdaHCfg6nvo9xpjRlscxSd0DF15YyqQBHwmIqtFZKbVgfyh/sTv0WUuNgFvhnMxrzcYGA9cXP9vu7q+56pCkPbQlVIqTGgPXSmlwoQWdKWUChNa0JVSKkxoQVdKqTChBV0ppcKEFnSllAoTWtCVUipM/H9eI9kt0jI/rAAAAABJRU5ErkJggg==",
	"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
	}