Slice sampling

slicesample.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 = 5.0\n",
    "    mu2 = -5.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 SliceSampling(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.w = 10.0\n",
    "\n",
    "    def sample(self):\n",
    "        u = self.inv_dist_fun(self.x) * self.rng.random()\n",
    "        x_min = self.x\n",
    "        while self.inv_dist_fun(x_min) > u:\n",
    "            x_min -= self.w\n",
    "        x_max = self.x\n",
    "        while self.inv_dist_fun(x_max) > u:\n",
    "            x_max += self.w\n",
    "\n",
    "        while True:\n",
    "            x1 = self.rng.uniform(x_min, x_max)\n",
    "            if self.inv_dist_fun(x1) > u:\n",
    "                self.x = x1\n",
    "                break\n",
    "            else:\n",
    "                if x1 < self.x:\n",
    "                    x_min = x1\n",
    "                else:\n",
    "                    x_max = x1\n",
    "\n",
    "        return self.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler = SliceSampling(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": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFpCAYAAACvXECGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/6klEQVR4nO3deXxU1f3/8deZZJLJTgIBQgIEEGR1QUAUUetWFZVqtdV+7W6p9au2fmsrdm9trf1+ta391aXa2mrrVvcFrBvuyL7IJgrIko0kQPZtkjm/P5IgO0mYmTN38n4+Hj4embk397wvJvnMOffec4y1FhEREXHH5zqAiIhIb6diLCIi4piKsYiIiGMqxiIiIo6pGIuIiDimYiwiIuJYl4qxMeZcY8x6Y8wGY8zsA2w/3RhTbYxZ0fHfz8IfVUREJD4lHm4HY0wCcBdwNlAELDbGPG+tXbvPru9Yay+IQEYREZG41pWe8RRgg7V2k7W2BXgMmBnZWCIiIr1HV4pxPrBtj9dFHe/t6yRjzEpjzEvGmHFhSSciItILHHaYGjAHeG/fOTSXAUOttXXGmPOBZ4GR+x3ImFnALIC0tLQTRo8e3b20IiIiHrZ06dJKa23uvu93pRgXAYP3eF0AlOy5g7W2Zo+v5xpj7jbG9LPWVu6z333AfQCTJk2yS5Ys6cYpiIiIeJsxZsuB3u/KMPViYKQxZpgxJgm4HHh+n4MPNMaYjq+ndBx3x5FFFhER6R0O2zO21rYaY64FXgYSgAestWuMMVd3bL8XuBT4jjGmFWgELrdaDkpERKRLjKuaqWFqERHpbYwxS621k/Z9vyvXjEVEpBcLBoMUFRXR1NTkOopnBAIBCgoK8Pv9XdpfxVhERA6pqKiIjIwMCgsL6bg9SA7BWsuOHTsoKipi2LBhXfoezU0tIiKH1NTURN++fVWIu8gYQ9++fbs1kqBiLCIih6VC3D3d/fdSMRYRkZhWVVXF3XffHfF2nn32Wdau3XfZhehQMRYRkZjW3WJsrSUUCnW7HRVjERGRg5g9ezYbN27kuOOO44YbbuDMM89k4sSJTJgwgeeeew6AzZs3M2bMGK655homTpzItm3buOWWWxg9ejRnn302V1xxBbfffjsAGzdu5Nxzz+WEE05g+vTpfPjhh8yfP5/nn3+eH/zgBxx33HFs3Lgxqueou6lFRKRbtnz5K/u9l3HeueR86UuEGhvZNuvb+23Puvhi+lxyMa27dlF8/Xf32jb0nw8dsr3bbruN1atXs2LFClpbW2loaCAzM5PKykqmTp3KRRddBMD69ev5+9//zt13382SJUt46qmnWL58Oa2trUycOJETTjgBgFmzZnHvvfcycuRIFi5cyDXXXMO8efO46KKLuOCCC7j00kt7+k/TYyrGIiLiGdZafvSjH/H222/j8/koLi5m+/btAAwdOpSpU6cC8O677zJz5kxSUlIAuPDCCwGoq6tj/vz5XHbZZbuP2dzcHOWz2J+KsYiIdMuherK+lJRDbk/Mzj5sT/hQHn74YSoqKli6dCl+v5/CwsLdjxClpaXt3u9gs0uGQiH69OnDihUrepwhEnTNWEREYlpGRga1tbUAVFdX079/f/x+P2+88QZbthxwESROOeUUXnjhBZqamqirq2POnDkAZGZmMmzYMJ544gmgvWivXLlyv3aiTcVYRERiWt++fZk2bRrjx49nxYoVLFmyhEmTJvHwww8zevToA37P5MmTueiiizj22GO55JJLmDRpEllZWUB77/pvf/sbxx57LOPGjdt9E9jll1/O//3f/3H88cdH/QYuLRQhIiKHtG7dOsaMGeM6RrfV1dWRnp5OQ0MDp556Kvfddx8TJ06MWvsH+nfTQhEiItKrzJo1i7Vr19LU1MRXv/rVqBbi7lIxFhGRuPTII4+4jtBlumYsIiLimIqxiIiIYyrGIiIijqkYi4iIOKZiLCIisofTTz+daD96q7upRUSkWwpnzwnr8TbfNiOsx/Mi9YxFRCTm1dfXM2PGDI499ljGjx/P448/zq9+9SsmT57M+PHjmTVr1u75qE8//XRuuOEGTj31VMaMGcPixYu55JJLGDlyJD/5yU+A9iUXR48ezVe/+lWOOeYYLr30UhoaGvZr95VXXuGkk05i4sSJXHbZZdTV1QHtyzqOHTuWY445hhtvvPGIz0/FWEREYt5//vMfBg0axMqVK1m9ejXnnnsu1157LYsXL2b16tU0Njby4osv7t4/KSmJt99+m6uvvpqZM2dy1113sXr1av7xj3+wY8cOoH3JxVmzZvHBBx+QmZnJ3XffvVeblZWV/PrXv+a1115j2bJlTJo0id///vfs3LmTZ555hjVr1vDBBx/sLvBHQsVYRERi3oQJE3jttde46aabeOedd8jKyuKNN97gxBNPZMKECcybN481a9bs3r9zjeMJEyYwbtw48vLySE5OZvjw4Wzbtg2AwYMHM23aNACuvPJK3n333b3aXLBgAWvXrmXatGkcd9xxPPjgg2zZsoXMzEwCgQBXXXUVTz/9NKmpqUd8frpmLCIiMW/UqFEsXbqUuXPncvPNN3POOedw1113sWTJEgYPHswvfvGL3UspAiQnJwPg8/l2f935urW1FQBjzF5t7PvaWsvZZ5/No48+ul+eRYsW8frrr/PYY4/x5z//mXnz5h3R+alnLCIiMa+kpITU1FSuvPJKbrzxRpYtWwZAv379qKur48knn+z2Mbdu3cr7778PwKOPPsopp5yy1/apU6fy3nvvsWHDBgAaGhr46KOPqKuro7q6mvPPP58//vGPYVkbWT1jERGJeatWreIHP/gBPp8Pv9/PPffcw7PPPsuECRMoLCxk8uTJ3T7mmDFjePDBB/n2t7/NyJEj+c53vrPX9tzcXP7xj39wxRVX0NzcDMCvf/1rMjIymDlzJk1NTVhr+cMf/nDE56clFEVE5JC8uoTioWzevJkLLriA1atXR6yN7iyhqGFqERERx1SMRUSk1yksLIxor7i7VIxFREQcUzEWEZHDcnV/kVd1999LxVhERA4pEAiwY8cOFeQustayY8cOAoFAl79HjzaJiMghFRQUUFRUREVFhesonhEIBCgoKOjy/irGIiJySH6/n2HDhrmOEdc0TC0iIuKYirGIiIhjKsYiIiKOqRiLiIg4pmIsIiLimIqxiIiIYyrGIiIijqkYi4iIOKZiLCIi4piKsYiIiGMqxiIiIo6pGIuIiDimYiwiIuKYirGIiIhjKsYiIiKOqRiLiIg4pmIsIiLimIqxiIiIYyrGIiIijqkYi4iIOKZiLCIi4liXirEx5lxjzHpjzAZjzOxD7DfZGNNmjLk0fBFFRETi22GLsTEmAbgLOA8YC1xhjBl7kP1+B7wc7pAiIiLxLLEL+0wBNlhrNwEYYx4DZgJr99nvOuApYHJYE0pYFc6ec8jtm2+bEaUkIiLSqSvD1PnAtj1eF3W8t5sxJh+4GLg3fNFERER6h64UY3OA9+w+r/8I3GStbTvkgYyZZYxZYoxZUlFR0cWIIiIi8a0rw9RFwOA9XhcAJfvsMwl4zBgD0A843xjTaq19ds+drLX3AfcBTJo0ad+CLiIi0it1pRgvBkYaY4YBxcDlwJf23MFaO6zza2PMP4AX9y3EEjtG79zCtJIPWJdTyPxBE1zHERHp9Q47TG2tbQWupf0u6XXAv621a4wxVxtjro50QAm/tGAjF298h58uepCvrn3JdRwRkV6vKz1jrLVzgbn7vHfAm7WstV878lgSSUsHjObiC27lmg+e5vKPXmdlvxGs6D/KdSwRkV5LM3D1IqW/+AVfXP86AMGERO455mJK0vpy1ZoXweoSvoiIKyrGvURLUTFV/36C9GDDp+8l+HlozLksHDAWf+iQN8KLiEgEdWmYWryv+plnwFqeH37KXu+/VXC8o0QiItJJPeNewFpL9XPPkXbSVCpSs/fbnhBqY0rZWlKDjQ7SiYiIinEv0PzRxwSLisg477wDbh9ZVcQvFzzAlLJ1UU4mIiKgYtwrGL+fPpddSvpppx1w+/rswVQnpXFC+fooJxMREdA1414hefgw8m655aDbrfGxtP/RTNr+ITYUwvj0GU1EJJr0VzfOtdXV07RuHfYwjy4tGXA0fVrqaVqroWoRkWhTMY5zDQsX8MnFl9C4fMUh91vW/2gA6ufPj0IqERHZk4pxnGtYugyTlERg/LhD7lednM6sM39A329+I0rJRESkk4pxnGtctozA+PH4kpIOu++2jAGYhIQopBIRkT2pGMexUFMTjWvWkHrCxC7tn1dfSekvf0nL5s2RDSYiIntRMY5jTatWQTBIysSuFWOftVQ9+hj1ixZFOJmIiOxJxTiOJY8ezeD77yN18uQu7V+c1o+E7Gwaly2PcDIREdmTnjOOYwkZGaRPn971bzCGwDETaFqzJnKhRERkP+oZx7GdjzxC07ruPTccGDuW5k2bCDU1RSiViIjsSz3jOFI4e87ur1ODjTw15xb+MeY8Hj/6zC4fIzB2LIkD+tNaVkZSYWEEUoqIyL5UjOPUiOoSADb0ye/W92WcdRaZZ58diUgiInIQGqaOU0dVFQHdL8bGmEjEERGRQ1AxjlPDq0uoDGRRnZzR7e+t+PNdFF13XQRSiYjIgagYx6khtdvZnDmwR98bqquj7u13sG1tYU4lIiIHomvGcer7p15LWrCx299XOHsO52xp4IbmZk7+7r8oTe+31/bNt80IV0QREemgnnGcavUl9miIGmBrxgCgvXctIiKRp2Ich8bs2My3P3iWrOa6Hn1/ZzEeWlsWzlgiInIQKsZx6JjKjXxu07u0+Hp2FaLBH2B+3niqktLDnExERA5E14zj0NDaMran9KHRH+jxMW458WvhCyQiIoeknnEcGlpTxtYe3km9J2NDYG0YEomIyKGoGMcZX6iNwXXlbOm47ttTJ5es4ukXf8zAhp1hSiYiIgejYhxnspvrqE5KY0vGkfWMdwYyCbQFGVqjm7hERCJN14zjzI6ULL587s+OeHh5a0Z/oP3xpoV548IRTUREDkI943h1hHNMN/hTqAxkUlBXHqZAIiJyMCrGceby9a9xw7LHw3KskrR+5NdVhuVYIiJycBqmjjPHVnxMUqg1LMeaN/iEHk2pKSIi3aNiHGcK6ipZkXtUWI71cuGJYTmOiIgcmoap40hyawv9mqopTs8N2zEzWuoJtDaH7XgiIrI/FeM4klfffn23JK3fYfbsmiE1Zfx77s+ZWromLMcTEZEDUzGOI0mhVj7MHrL7saQjVZrWlxCGQfW6iUtEJJJ0zTiOfJQ9hBtOuz5sxwsm+KlI6aM7qkVEIkw9Yzmk4vR+5NdVuI4hIhLXVIzjyC/e/xvXL38irMcsSe/XPkytBSNERCJGw9RxZGRVEYsHjAnrMd/MP54NWQX4sIQ4slm9RETkwFSM40RbXR05zbUUp4fnTupOa/oNZ02/4WE9poiI7E3D1HGiZcsWIHyPNXXyhdoYXlVMv4aqsB5XREQ+pWIcJ4LbtgHtjyOFkz/Uxl1v/oGzti0J63FFRORTKsZxIiEnh3cGHUNZmItxc2ISO5MzGNCwM6zHFRGRT+macZxImzKFW6d8JSLHLk3ry8D6HRE5toiIqGccN2xLS8SOvT01h4HqGYuIRIyKcZzYNPNzYVvHeF9laTnkNlSREGqLyPFFRHo7DVPHARsKESwqonro0Igc/42CiazNKYzIsUVERMU4LrRWVGCDQbanZkfk+EUZ/SkK0+ITIiKyPxXjOBAsLgZge1pORI7vC7UxefuHbE+NzPFFRHo7XTOOA8GiIgDKUsP7WNNuxvDjRQ/xmaJlkTm+iEgvp2IcB5KGDSPnq1+lPELD1CHjozw1W88ai4hEiIpxHEiZMIEBN8+mJcEfsTbK0nLI07PGIiIRoWIcB4IlJYSamiLaRllqjnrGIiIRohu44sCWK79MygknQMJpEWujLK0vWS0NtNXWkpCREbF2RER6IxVjj7PBIMGyMjIL8qE0cu28PvgEFg0Yw1spKZFrRESkl+rSMLUx5lxjzHpjzAZjzOwDbJ9pjPnAGLPCGLPEGHNK+KPKgQS3b4dQiKSCgoi2syuQydbMgZhEfX4TEQm3wxZjY0wCcBdwHjAWuMIYM3af3V4HjrXWHgd8A/hrmHPKQXQ+1uTPj2wx9oXamLHpPRqW6fEmEZFw60rPeAqwwVq7yVrbAjwGzNxzB2ttnbXWdrxMAywSFbuLcYR7xiHj46o1L1L7yqsRbUdEpDfqSjHOB7bt8bqo4729GGMuNsZ8CMyhvXe8H2PMrI5h7CUVFRU9ySv7SDnuOPrPvgn/wAGRbcgYKlL6ECyN4IVpEZFeqivF2Bzgvf16vtbaZ6y1o4HPAbcc6EDW2vustZOstZNyc3O7FVQOLPmoo+j7ta9F5VpueWo2wZKSiLcjItLbdKUYFwGD93hdABz0L7K19m1ghDGm3xFmky5oXLUqagWyPEXFWEQkErpSjBcDI40xw4wxScDlwPN77mCMOcoYYzq+nggkAZquKQqKrrueijv/FJW2ylOzaduxI+ITjIiI9DaHLcbW2lbgWuBlYB3wb2vtGmPM1caYqzt2+zyw2hizgvY7r7+4xw1dEiE2GKS1vJzEQXlRae/54dMYtWQxvkAgKu2JiPQWXbrQaK2dC8zd57179/j6d8DvwhtNDqe1vBxCIfx50SnGDf4UEtLTo9KWiEhvormpPazzzmZ/3qCotJfc2kz5HXdQv2BhVNoTEektVIw9bHcxzo9OMW71JbLjbw9Qv3BBVNoTEektVIw9LHXKFPLvvBN//n6PfUdEmy+BxAEDaNUd1SIiYaWJhj3MP2AA/s+eE902Bw0iWKxiLCISTirGHlI4e85er48v/4h6f4CPsodELYN/0CAaNT+1iEhYaZjaw761+nm+uP71qLbpHzSIUF0dNhSKarsiIvFMxdjD+jdUUZGaHdU2c6+7lpEL3sf49KMjIhIuGqb2qNRgI2mtTZSn9Ilqu1rPWEQk/NS98aj+DVVA+xSV0dS6axfFP/gh9fPnR7VdEZF4pmLsUbmNuwCoiHLP2JecTM0LL9C4anVU2xURiWcac/So1X2H893TrmdLRoTXMd6HLzWVhGyt3iQiEk4qxh7VGOVHmvbkz8tTMRYRCSMNU3vUySWrOLlklZO2/fmDVIxFRMJIPWOPumTDW7QZH/MHTYh620kjRhDcXh71dkVE4pV6xh6V21gV9TupO/X/3vcY9u/HnbQtIhKPVIw9yBdqo29TDeUpboqxiIiEl4qxB+U01ZJgQ1Sk9nHSfktRMVuu/LKeNRYRCRMVYw/q3/GMsauesS81hYYlS2jeuMlJ+yIi8UY3cHnQupyhfOncn9GQGHDSfkJ2NiY5mWBpqZP2RUTijYqxB1njY1cg01n7xpj2Z41L9XiTiEg4aJjag07ftozPbXjbaQb/oDxaS9QzFhEJBxVjDzq9aDlnblvqNEPKccfhH+pmBjARkXijYWoP6t9YRWlqjtMMuddf77R9EZF4op6xB+U27KLC0YQfIiISfirGHpMabCS9tSnqSyfuq3HNGjacdTYNixc7zSEiEg9UjD0mp6mWVuOj3NGEH50S0tMJFhXRUlzsNIeISDzQNWOPKcroz8yLbsNY6zRH4sCBALTqWWMRkSOmYuxBIeMD4zaDLzmZhH79COrxJhGRI6Zi7DFnbVnMyKpt3HPsJU7aL5w9Z/fXf2xNoe6tlfxkj/c23zbDRSwREU9TMfaYiRUfMXrnFu5xHQRYkDeOhFCb6xgiIp6nYuwxuY1VVDq+k7rTY0ef5TqCiEhc0N3UHpPbsIvtMfSMsbEhjA25jiEi4mkqxh7isyH6NdU4f8a40wnbP+S5F25meLVu4hIRORIqxh6S3tJIZSCLMsdTYXaqTk7HH2ojt2N9ZRER6RldM/aQmuQ0vvbZH7uOsVtnD71/Q5XTHCIiXqeesfRYdVIazb5E+qtnLCJyRFSMPeSsLYv51fz7SQy1uo7SzhgqUrPJbVAxFhE5Ehqm9pCRVUWM2bmFVl/s/G97aeiJ1CWluo4hIuJpsfNXXQ4rt7Eq5pZOfHrk6a4jiIh4noapPSS3sYryGHmsqZOxIbKbavDpWWMRkR5TMfaQ3MZdVDheOnFfZ25dyiP/+RX9dd1YRKTHVIw9wra2si1jAJ9k5rmOspfODwe6iUtEpOd0zdgjTGIiP5j+365j7Kc8pf0a9oDGXaxynEVExKvUM5YjUpmSBUCuJv4QEekxFWOPqH7hBe6edztZzXWuo+wlmOBnZ3IGuY1VrqOIiHiWhqk9ouWTzQyt2U6dP8V1lP38a8xnY2olKRERr1Ex9ohgaSk7A5m0+RJcR9nPS4VTXUcQEfE0DVN7RLCslPIYe6ypU2qwkVG7toK1rqOIiHiSirFHtJaUUpESm0PBZ29dwp1v/YnMlgbXUUREPEnD1B6RcsIJrNroOsWBdS6lqHWNRUR6Rj1jjxh062+YM+xk1zEOqLzj5i2taywi0jMqxh5gY/xarHrGIiJHRsXYA2pfe431k6cwuHa76ygHVJ2URlOCn/561lhEpEdUjD2gtbSUUG0t1UnprqMcmDHcMfFyXht8guskIiKepBu4PCBYUooJBKhJSnUd5aDezT/WdQQREc9Sz9gDgmVl+PPywBjXUQ5qQP0OppStdR1DRMSTulSMjTHnGmPWG2M2GGNmH2D7fxljPuj4b74xRt2kMAqWlrQX4xh25ral/HLBA4RaWlxHERHxnMMOUxtjEoC7gLOBImCxMeZ5a+2e3aBPgNOstbuMMecB9wEnRiJwb5Rxxpkk9OkDy10nObjOpRRby8pIGjLEcRoREW/pSs94CrDBWrvJWtsCPAbM3HMHa+18a23ncy0LgILwxuzd+n17Ftlf/ILrGIfU+axxsKTUcRIREe/pSjHOB7bt8bqo472D+Sbw0oE2GGNmGWOWGGOWVFRUdD1lL2ZbWmirq3cd47AqOubNDpaqGIuIdFdXivGB7ho64CwUxpjP0F6MbzrQdmvtfdbaSdbaSbm5uV1P2Ys1rFjBR5MmUb9goesoh1QZ6AO0X98WEZHu6UoxLgIG7/G6ANjvL64x5hjgr8BMa+2O8MST1o6epn/gAMdJDi2YkMhN066mz+cvdR1FRMRzulKMFwMjjTHDjDFJwOXA83vuYIwZAjwNfNla+1H4Y/ZencO+iQMHOk5yeB/kHoV/QH/XMUREPOewxdha2wpcC7wMrAP+ba1dY4y52hhzdcduPwP6AncbY1YYY5ZELHEvEywtIyEnB18g4DrKYY3atZWqZ551HUNExHO6NAOXtXYuMHef9+7d4+urgKvCG03AG88Ydzq1aAVlv7yPrM/NxMTwBCUiIrFG02HGuKyLZmKDQdcxuqQiNRvb1ERbVRWJ2dmu44iIeIaKcYzLumCG6whd1jnxR7CkRMVYRKQbNDd1DAs1N9O8YQOh5mbXUbqkvONZ41Y9aywi0i0qxjGs+eMNbLrgQurfecd1lC6pSNEsXCIiPaFh6hgWLC4GwJ9/qAnPYkdNUirDX3zBM3lFRGKFinEM212MBw1ynKSLjCH5qKNcpxAR8RwNU8ewYEkJvvR0fJmZrqN0We28eex86CHXMUREPEXFOIYFi4vx5+d76pndujfeoPL++13HEBHxFA1Tx7Ccr32NUH3sr9i0p8S8PNoqKgm1tOBLSnIdR0TEE1SMY1jaiVNcR+g2f1779e3WsjKShgxxnEZExBs0TB2jQo2N1L3zLq27drmO0i2dU3fq8SYRka5TMY5RzZs2se1b36Jh8WLXUbrFP6i9GLeWb3ecRETEOzRMHaOCJe1LRnvtmV1/QQGjliwhIT3NdRQREc9QzzhGdT5jnOSxYmx8PhVikRhjW1poXLPGdQw5BBXjGBUsLsGXloYvK8t1lG7b+fDDVN577+F3FJGIq339dT4+80xKbvzB7vdqXn3Vc09qxDsV4xjlxWeMOzUsXkL1c8+7jiHS61U9/QxF116HP7c/A350M9D+t6X4hv9h6zevoq1OBTlWqBjHqP7f/x8G/vIXrmP0iD8vj2BpKdZa11FEeq36RYso/clPSDvpJIY+8jDp06cD7feh5P/h9zSuWkXpT3+i39MYoWIco5JHjCD1+ONdx+gRf14etqmJNo89liUST6qffoakwYMp+H9/whcI7LUt8+yzyf3ud6l96T9UPfmko4SyJ91NHUMKZ88BILm1mdOLVrA8dyTlaTmOU3WfP7994o9gcTGJOd7LLxIP8m79Da0VlfjSDnxDZd+rvkn9u+9SfvsdZH72syR4aA78eKSecQwaVL+D7614glFV21xH6RF/QQG+jAzaqqpdRxHpdYLFxQTLyzE+H/4B/Q+6n/H5GPjTn5D3q1/hy8iIYkI5EPWMY9CAhp0AbE/1Zq8yedQojl68yHUMkV5p++2307h0GedO/T5tvoSD7rf5thkkjxxJ8siRUUwnB6OecQwa0NB+rdWrxdiLd4CLxIPmTz6h9j8vkzVz5iEL8b4q77mH8jvvjGAyORwV4xg0oGEnjQlJ1CSluo7SYxX/78+U/eZW1zFEepWdD/wdk5REzle/0q3vaykqYucDf/fcXPjxRMU4Bg1o2EV5ajZ4uIfZsnkzdW++6TqGSK/RVlND9QsvkHXRhST269et7835ylewzc1UP/tchNLJ4agYx6A/HP8FfnXi11zHOCL+wQUES0qwra2uo4j0Cg2LFmGbm+nzxcu7/b2Bo48m5fjjqXrsMT137IiKcQyqS0qlJD3XdYwjklRQAG1tBMvKXEcR6RUyzjqLo96YR8r4cT36/uwrLqdlyxYaFiwIczLpChXjGJPe0sCX1/2HwbXeXoLQXzAYgOA2bz6eJeIlnb1Z/8CBPT5Gxmc/S9Yll5CguQGcUDGOMYNry/nS+tcYWL/TdZQjkjRkMMmjRmFb21xHEYl723/7W4qu/+4RDTH7kpMZdOtvCBx9dBiTSVfpOeMYk1e/A4DStL6Ok/RM5yxiAIydBXOqYc6n722+bYaDVCLxy7a2UvPiHFKnTAnLY4VNH32EbQn2eLhbekY94xiTV19JCOPZZ4xFJLrqFyykbedOMmecf8THstZS9J1rqPjjH488mHSLinGMyavfQWVKFsEE7w9aXLX6BX7x/t9cxxCJazVz5+JLTyf91FOP+FjGGDIvuID6+fNprawMQzrpKhXjGNOvqdqzQ9T7CrS2MHrXFtcxROJWqKWF2ldfJeOss/AlJ4flmJkzzodQiNrXXg/L8aRrvN/9ijOzp11NSmuz6xhhUZqWQ1ZLA6nBRhr8Ka7jiMSfYJC+V11F6uRJYTtk8siR+IcOofa118i+/IthO64cmnrGscYYGv2Bw+/nAZ3Xvb1+Z7hIrPKlpdHv27NInTgxbMc0xpBx1lk0rlhBqDk+OgZeoGIcQ/LqK7lh2eMMqYmPiTI6h9sHNqgYi4SbbWuj5tVXaaurD/ux+151FSPfejNsQ99yeCrGMaSwuoxzti4muS3oOkpYlKX2ZdGA0TTESU9fJJY0fvABxdddH5E54BOzs/GlpYX9uHJwumYcQ/Lq2+9ejJcbuOqTUvj5SVe5jiESl+reeBMSE0mffkqPvn+vOQEOYPWMTHb85T4G/+2v6iFHgYpxDMlr2EGtP4U6Dy+deCA+GyJkNAgjEk5LH32Oqj5DOe+370asjYYlS2hYsID0006LWBvSTn8hY0he/Y646RV3+vYHz3H/a79zHUMkrrRs20Zh7XYWDozcLFmpU6ZgUlKoe+vtiLUhn1IxjiHGwrb0/q5jhFWdP8DA+p344+Q6uEgsqJ//PgAL88ZGrA1fcjJpJ55I3TvvaFnFKNAwdQz58bRZriOEXUl6Lj4sefU72JrZ8xVlRORTfb5wGRe+XU9pWr+ItpN26nTq3nyTls2bSR42LKJt9XbqGUtEFae3/7EYVK+p9UTCxRhDUUbkR9HSTz2NtFOnY5uaIt5Wb6diHCMaFi/m1vf+svuO6nhRnJYLQH5dheMkIvGhfuEiSmbfTJ+m2oi3lVSQz5D77iMwZkzE2+rtVIxjRNO6dRxf8TGNCfH1CEF9UgrPDp/OJ1mDXEcRiQu1r79GzUsvUR/F5/dbKys1G1eEqRjHiJbNm6lLDFCVnO46Stj95ZiZLOuvBctFwqH+vfmkTppEMMEflfYali/n41OmUz9/flTa661UjGNE8yefUJyeC2FYHDzW+GyIfg1VrmOIeF6wtJSWjRtJmzYtam0Gxo3DpKZS97YecYokFeMY0bJ5C0Xpua5jRMSlH7/BP1/5NYE4WY1KxJXO3mk0i7EvKYm0KVOof08940hSMY4Btq2N5GHDWJ8zxHWUiCjpePxiUF183ZwmEn2GlOOPJ3nUyKi2mjZtGsGtW2kpKopqu72JinEMMAkJDHngb7wwvGdzzMa63cU4zu4UF4m2Pp+/hMJHH8FE+XJW2rSTAdQ7jiAVY4m4ko5njfPVMxbpsVBTEzYUctJ20rBh5N16K+mnneqk/d5AxTgGVN5/PxtnXEBCqM11lIhoSkymMpBJQV256yginrXzHw/y8bRTCNWHf/3iwzHG0OeSi/EP1Cx6kaJiHAOaP/6YUEMDbb4E11Ei5sGx5/H6kEmuY4h4Vv38+SQOHOhsneG2ujqqnnyS5k8+cdJ+vFMxjgEtn2wmqXCo6xgR9dqQyazIje5NJyLxIlRfT8Py5aR3XLt1wTY3U/qTn1L78svOMsQzFWPHrLW9YhL25NZmxlVuoq2uznUUEc+pX7wYgsGoPtK0r8S+fUkeM0Y3cUWIirFjbTt2EKqtJamw0HWUiBq1axu3v3s3jStXuo4i4jn1783HBAKkTJzoNEf6tJNpWLHCyXXreKclFB2zra1kXXwxKcccA2uLXceJmG0ZAwBo2bgRHH66F/GKwtlzdn89bkcmQ46ewUs/f81hIkg7+WR2/PVv1C9eTMbppzvNEm+61DM2xpxrjFlvjNlgjJl9gO2jjTHvG2OajTE3hj9m/PIPHMig395KynHHuY4SUVXJ6dT4U2neuMl1FBHPWdN3GC8NO8l1DFJOOAETCND84XrXUeLOYXvGxpgE4C7gbKAIWGyMed5au3aP3XYC1wOfi0TIeLHnJ91OqcEmGhKT43JO6r0Yw7aM/gzYuMF1EhFPOaqqiOS2IGtzhmJN9K4sHujvFUD6GT+ibnMqm6OWpHfoyv/ZKcAGa+0ma20L8Bgwc88drLXl1trFQDACGePaLxY8wG/m3+86RlRszRhAi3rGIt1y6cdvMHvxP7HExgf2uqRU1xHiUleKcT6wbY/XRR3vyZGylqE1ZZSl5bhOEhXPjTiFgrvvwlrrOoqIJ/hsiOPLP2ZZ/1ExM3qW3tLAzxb8nZqXXnIdJa50pRgf6CegR39NjTGzjDFLjDFLKioqenKIuJLdXEtmsIEtHTc3xbstmXmkHn981OfVFfGqo6qKyAw2sDx3lOsou9X7A4zetYXaeW+4jhJXulKMi4DBe7wuAEp60pi19j5r7SRr7aTc3PhcLrA7htZuB2BrZu+YYs4XaqP6xTk0rlrlOoqIJ0ws/wiA5f1jZ8Ica3ysyB1J/fvva5QrjLpSjBcDI40xw4wxScDlwPORjdU7DKkpA+g1PeOQ8VH2s59R/fwLrqOIeMKEyk1syMqnOjnDdZS9LMsdRVtlJc0ffeQ6Stw47N3U1tpWY8y1wMtAAvCAtXaNMebqju33GmMGAkuATCBkjPkeMNZaWxO56N73Yc5Q/jn6HHbF2C9axBhD0lFH0fzxx66TiHjCz0/6BjmNsfdndHn/9mHz+nffI3D00Y7TxIcuTfphrZ0LzN3nvXv3+LqM9uFr6YaPsofwUfYQ1zGiKnD0KGpffQ1rra4dixxGqy+R8hi8wXNHShaZ559PQt/Yy+ZVmg7TFWsZvXMLgdZm10miKvno0bRVVdFaruUURQ7l8x+/wZc+fMV1jIPK//0d9Pnc51zHiBsqxo7kNNXwh7f/H2dtXeI6SlQFRrcPaTWv1ww+Iody3uaFjNq17fA7OhRqaqJ11y7XMeKCirEjw2vab0jf0kvupO4UmDCBo15/jbTp011HEYlZLdu2kV9fybL+sXs91gaDfHza6ez4y32uo8QFFWNHRlS1F+ONWYMcJ4kuX3Iy/vx8XS8WOYT6994DYGn/2Hm+eF/G7ydl/Hjq3n3HdZS4oGLsyPDqYkpT+9LgT3EdJepq571B+e23u44hErPq33uP7SnZFKfH9nwMadNPoWXDRoIlPZp6QvagYuzIiOqSXtcr7tS0ejU7Hvg7oaYm11FEYpIvM5N38o+JmSkwDya943JT3bvvOk7ifSrGjvx+4hd5YtRnXMdwInn00RAK6XljkYMY9Jvf8LfxF7qOcVhJw4eTmJdH/TsqxkeqS88ZS/it7TvMdQRnAqNHA9D04YekTJjgOI1IbAk1NOBL9cbKSMYYBv7spyT27+86iuepGDsweucW+jZWM3/Q+KiuTxor/AUF+FJTtUC5yAFs+frXScovgKQzXEfpkozP9M4RvnDrfZUgBnx2y0KuW/lUzKxPGm3G5yMwdiyhulrXUURiSltVFU2rVpM0zFsjZ3VvvUXtvHmuY3iaesYOHFVVzMas/Ji/OSOShjz4D0xCgusYIjGl7t33IBQiffop8Fix6zhdVnn//diGRjLO8EZvPhapZxxlSW1BCmtK+Sh78OF3jmMqxCL7q3vzTRJycgh47F6K9FOm07R2La2Vla6jeJaKcZSNqCom0YZY38uLcVtNDVu/8U2qX3jRdRSRmGBbW6l/5x3Sp0/33IfVtOmnAJ9OViLdp2IcZSOrigBY38tWa9qXLyODxjVraFi0yHUUkdgQCjHgRzfT5/Ivuk7SbYExY0jo25e6tzUbV0/pmnGUPT98GosGjmFXINN1FKeMMaSMG0fj6tWuo4jEBJOURNbMma5j9Ijx+Ug/7TQaP1ip5VF7SMU42oyhLK2v6xQxITBhAjv++ldCTU34AgHXcUScqnrmWdKmTMafn+86SpcUzp6z1+vU4EQax03F3jx393ubb5sR7ViepWHqKMpoqeeGZY8zvMo7d0lGUsqE8dDWRtPada6jiDjVUlRM6c03U/Pqq66j9FiDP9Ar500IF/WMo2j0zq2cs3Uxrw2Z5DqKM3t+ms5uquGH/Ubw/XveY23f9onm9UlaeqO6t94EIOP0053mOFLnfzKfM7Yt48bp/92rH93sCRXjKBq9cwttxsfHfQpcR4kJuwKZ3HzKd1zHEHGu7vXXSSosJKmw0HWUI9JmEhi3czOFNaVs7qUL4fSUxhSiaMKOTWzIyqcpMdl1lJgSaG0Ga13HEHGiddcu6hcuIuPss11HOWKLBo4hhOGk0jWuo3iOesZREmpp4ehdW3lh2Mmuo8SUk0tW8aPF/+TqM26kKEOTzUvv0XnJZnLZOn7eFuK/Pkrl431uivKaXYFM1mcPYWrZGh4d7f0PF9GknnGUtJaVsT01m9X9hruOElO2ZA4kwYYYt+MT11FEnFg8cAz/dd7P4+by1YK8sYyqKiK3YZfrKJ6innGUJA0ZwqyzbtJw7D6K0/qxKzmd8Ts28XLhia7jiDhRnZzuOkLYvDPoWAbW78RnQ66jeIp6xlFiO4uw7jDcmzGsyRmmnrH0StOLV/Dbd+8lq7nOdZSwKU3vx5+Ov4ztmk+hW1SMo8C2trLxrLM5/5P5rqPEpFX9hpPXsJP+9TtdRxGJqlOLVlJQV05NUqrrKOFlLSN3baO1osJ1Es/QMHUUNH6wimBxMbWD4uwXLkwWDhxLgg3pLnPpVQKtzUwq/5BXhkyOu8ky+jVW86e37qTqmQT6zfqW6zieEF8/ATGq/r33wOdjee4o11Fi0va0vjxz1GnUJKe5jiISNSeXrCbQFuStguNdRwm7ytQ+rMseSs1LL7mO4hkqxlFQ/+67BCaMpy7ehqLCKL2lgenFK7Btba6jiETFGUVLKUvNZm1OoesoEfF2/rE0r1tH8ybdD9IVKsYR1lZdTeOqVaRPm+Y6SkybvP1DfrT4X5qnWnqNRQPG8O9RZ8btTZ3v5B8LxlAzx9vPTkeLinGE2WCQvt/4OhlnneU6SkxbnjsS0OLk0ns8P2I6LxVOdR0jYnakZJE69URq581zHcUTVIwjLLFfP/rfeCOBsWNdR4lpVYEMPs7Kp+7NN11HEYm4mldfJSXY5DpGxOXd8msKH33EdQxPUDGOINvWRv2CBdhg0HUUT3g/bzyNK1fqcQiJa01r11J83fWctW2J6ygRl1SQr7XKu0jFOIIaV6xg69e+Tu1rr7mO4gnv540Ha6lfuMh1FJGI2fXY45hAgHkFJ7iOEhX18+ez6eJLaKuudh0lpuk54zAq3GeS92+tep4LfQlMfa2Bhrd0E8PhbM4cyIiX/0PS0KGuo4hERFtdHdUvvkjmjPOpb0txHScqfJlZNK9bR/WcOeR86Uuu48Qs9YwjxVpOLl3FstxRNPg1TNMlxqgQS1yreeEFbEMD2Zdf7jpK1ATGjSUwdiy7Hnnk02mBZT8qxhEyorqEgQ27mJ833nUUT2mrraX4f/5HkwVIXGpcsYLA2LEExveevwvGGLK/8mVaNmykfr6mBD4YDVNHyMklqwiaBBaoGHeLLy2NxpUf0FZdQ+Z557mOI9Jj+162AsCcSlrhZOpvnhv9QA5lnn8+5bffwc6HHtKcCwehYhwhD485h3fzj9EUj91kfD4yL7qQHX+5j2B5Of7+/V1HEgmL1GAjDf4U6pN6x7XiPfmSkhjwgxsxKb3v3LtKxThCQsbHJ1mDXMfwpKyLLmLHPfdS8+Ic+n7j667jiByx4VXF3PHOn7llytdYNuBo13GiZu/RgUQgCPM+fW/zbTOinilWqRhHwDUrn6YyJat9qjvpls5f3j9kD+aTex/imvW5e00XqF9e8aJLN7xJyPhYnz3EdRSn0loauWTDW/yn8EQqUrNdx4kpuoErzHIaqzl380KymutdR/G0Z0acxsKB4/CHtHCEeNugugpOLVrB3MKpvXKIek+prU1c9vEbfPEjTZG5LxXjMLvwk/kk2BAvDj/ZdRRPe7vgOB4cex7BBA3eiLd9Zd3LBBMSefqo01xHca4iNZtXhk7hnC2L6N+w03WcmKJiHEbJrS2ct/l9FuSNozStn+s4nuezIU4sXUO/hirXUUR6pF9DFaeUfMDTI05jVyDTdZyY8NioM8HAlz7UzIR7UjEOozO2LSWrpYFnRkx3HSUu5DTW8NNFD/L5DW+6jiLSI5WpfbjmM//DUyPVK+5UmdqH54dN4+yti2lcvcZ1nJihYhxG63KG8tRRp7G673DXUeJCZWofXh0yifM3L6BfY5XrOCLdkt7SAMDWzIE0+Hv3teJ9PTL6bN7JPxZfWqrrKDFDxTiMNmcN4q/jL4zbxcJdeHTUWfhsiC/ohg/xkLa6Ou6Zdztf+vAV11FiUoM/hdsmX0nysGGuo8QMFeMwaKurp/QXv6B/vW5ICLfytBxeHjqFczcvJL9OSyuKN5T/7n/JaaplyYDRrqPEtGB5OcXfv5FgWZnrKM6pGIdBxZ/upOrxf9Onpc51lLj08OhzKE3vS06jlmCT2Fc7bx5VTzzBkyNP56Ne/lzx4dimJmrnzaP0xz/BhkKu4zilYnyEGj/4gF3//BfZV1yuX7wI2RXI5OozbmRV7lGuo4gcUrCkhNIf/4TkMWP455jPuo4T85KGDGHATT+k/r332HHffa7jOKVifARCDQ2U/vjHJObmknvDDa7jxDVrfPhsiMp7/0KwpMR1HJEDatm8GZOURP4dd9Dq0zPyXdHni18k84ILqLjzT9S9847rOM6oGB+Biv/3Z5o3bCTv1ltJyMhwHSfu9WusYsf991P8/RuxwaDrOCK7da7Tm3byyYx45WWSh+vGpK4yxpD3q1+SPGoU5bff0WuHq42rxZ4nTZpklyxZ4qTtcGndtYv6+fPJmtE+X/IBl0yTsFo5HUq+fyNZl36evFtuwejOdXGk8/fdZ0PcuPRRVvcdxtxhmnmvp3Iaq7HG7DU5SjzORW+MWWqtnbTv++oZ90D9/PmEWlpIzM7eXYglOrJmzKDv1d+m+smnqPj9H3D1YVIE2mfd+8nCB/lM0XLSgk2u43jazpQsdgUy8dkQ1654kjE7NruOFFW6qNENhTe9yMUb3+aq1S/y2Kgz+efYc11H6pVyv/td2nbuYueDD9Ln0s+TNHSo60jSC+XVVzJ78b8YUVXMXcdczIvDp7mOFBcyWuo5rmID52xdzK7Hc+nzhct6xQiYhqm7qK2unj9f8h3O2bqYdwdN4PaJV9CcmOQ6Vq/TOWxlQyGaP/qIwOj25zjbqqtJyMpyGU16kba6OpaedCoAd0y8nIV54xwnii8ZLfX8cMkjTCpfz/yB47jn2IupTOmz335eHMY+2DC1inEXNK5aTfENN9BcVMzjo87gn2M+izUa4Y8Vp29bxnc+eJYHx57LXQ//EpOoAR8JP2stTavXkDJhPADXfuFHLBo4hvLUHMfJ4pPPhrh4w1tc+eErFKX357rTv7ff7IbxVIy7VFGMMecaY9YbYzYYY2YfYLsxxvypY/sHxpiJ4Qjtkm1p2T0rTGK/vvgyM/jh9Gt4aOx5KsQxZkOfArZkDuC6lU+z6YIL2fnww7TVaT1pCY+2qip2PfEEn3z+82y+7LLdixu8OHyaCnEEhYyPp0Z+hu+ccSN/HX8BGENKsImrVj3PsOr4e7zxsD1jY0wC8BFwNlAELAausNau3WOf84HrgPOBE4E7rbUnHuq4sdYzttYSLC6mccVK6t5+i/q33iYwbixDHnhg9/ZhN891nFIOylpOLl3NbxqX0bRqFSkTJ1L4yMMAtGzbhj8/H+PThyg5vM7f9f4NO7l+xZMcW7GBRBtia3p/nj7qNOYNnkgwwe86Zq90XMXH/Gr+X/HbNkrS+jJu5jmkTp1K2snTSEhPcx2vSw7WM+7KeN4UYIO1dlPHgR4DZgJr99hnJvCQba/sC4wxfYwxedba0jBk7zJrLQSD2GAQ29qKbW0lIScHYwytlZUES0tpq64hVFNNW00toYYG+n7j6wAUX389ta+2r6+Z0KcP6aefTuYFF+w+dm+4gcDTjGH+oAl8hgmMzt5CarCJZbPnkNQW5KkXf0xzQhLbMnKZMu1Y/IMLSJ8+ndSJE7HBIE0ffohJTsYXCGCSA/iSk/ClpWH8/vafqVAIfD79DERBZ+eg89/atrWBtZ/+174R428vhqHm5vb/P53brIWEBHyBQPv2+npsayuhpiZCDQ3YpiZ8qakkDR2KtZaqxx6jtaKC1opKgtvLaNm4qeP3fiQ1SWnkNNXwzIhTeavgODZm5WsRGMdW5I7kv877GacWr2Ry2TqyH3+KlEce5Rtnz6Y0rR+nFS3npNI1lKf2YUcgi1/+11QSMjNJP/VUTGIirTt3Yhsbwe/H+P0YfxK+JD8myf39P10pxvnAtj1eF9He+z3cPvlA1IpxxZ/+ROXd9+z3/tErlmMCASr/ch+7/vnPvTcmJpLz5Ssxfj+ZF17IreWZrM8ewidZgwgZH8yphjl6dthrPsz59O5qYy13HncZR1dtZVBdJY3Ll1Mzdy4JaWmkTpxIsKyMzZd9Yb9jDPjZT8n50pdo/vBDPrn4kk83+Hzg8zHo1t+QddFFNCxbztZvfnO/78+/4w4yzvgMde++R9F11+23ffA9d5M2dSo1L79Cyez9rvww9KEHSZkwgapnnqXsllv22z7siX+TPGIEOx9+mPI7fr/f9hFz5+AfOJDK+++n8p5799s+8s03SMjMpPzOO9n54EP7bT960UJMYiLbf/tbqp54cq9tJjmZUe/PB6D0pz+les7cvYphQk42I+e1r7JV9N3vUdvxdWdBTSooYMTL/wFg6ze+Qf389/c6fvLRRzP8uWcB2Hz5FTStWrXX9j1HPT65+BJaNm3aa3va9OkMub99asWNMy6gdZ9FCDLOO5dTks8C4IkX/5fU1maqktPZGchkW0Z/5i9rgHxoSkzmmjNu3O/fRtyqTUpjzrCTmTPsZPxtrQyrKaGs43JBVnM9R1UVcXLpKvyhNkpnPwfAjIt+R8iXwH+vfIoLPtn7580kJzN65QoASm66ieoX57R/+Aby77yTjDM+E5Xz6sow9WXAZ621V3W8/jIwxVp73R77zAF+a619t+P168APrbVL9znWLGBWx8ujgfXhOhGgH1AZxuO5pHOJTfFyLvFyHqBziVXxci6ROI+h1trcfd/sSs+4CBi8x+sCYN+r513ZB2vtfUBEZgM3xiw50Di8F+lcYlO8nEu8nAfoXGJVvJxLNM+jK3e0LAZGGmOGGWOSgMuB5/fZ53ngKx13VU8FqqN9vVhERMSrDtsztta2GmOuBV4GEoAHrLVrjDFXd2y/F5hL+53UG4AG4OuRiywiIhJfujQ7grV2Lu0Fd8/37t3jawv8d3ijdVs8LYapc4lN8XIu8XIeoHOJVfFyLlE7D2czcImIiEg7zYIgIiLiWFwVY2PMccaYBcaYFcaYJcaYKa4zHQljzHUd05CuMcb8r+s8R8oYc6Mxxhpj+rnO0hPGmP8zxnzYMeXrM8aYPq4zddfhprb1CmPMYGPMG8aYdR2/H991nelIGGMSjDHLjTEvus5yJDomfHqy4/dknTHmJNeZesoYc0PHz9ZqY8yjxphAJNuLq2IM/C/wS2vtccDPOl57kjHmM7TPbHaMtXYccLvjSEfEGDOY9ilVt7rOcgReBcZba4+hfYrYmx3n6ZaOqW3vAs4DxgJXGGPGuk3VY63A9621Y4CpwH97+FwAvguscx0iDO4E/mOtHQ0ci0fPyRiTD1wPTLLWjqf95uXLI9lmvBVjC2R2fJ3FAZ519pDvALdZa5sBrLXljvMcqT8AP6T9/5EnWWtfsda2drxcQPvz9F6ye2pba20L0Dm1redYa0uttcs6vq6l/Y9+vttUPWOMKQBmAH91neVIGGMygVOBvwFYa1ustVVOQx2ZRCDFGJMIpBLhehJvxfh7wP8ZY7bR3pP0VM9lH6OA6caYhcaYt4wxk10H6iljzEVAsbV2pessYfQN4CXXIbrpYNPWepoxphA4HljoOEpP/ZH2D6ohxzmO1HCgAvh7x5D7X40x3li9YR/W2mLaa8hW2qd1rrbWvhLJNj238Ksx5jVg4AE2/Rg4E7jBWvuUMeYLtH9COyua+brjMOeSCGTTPgQ3Gfi3MWa4jdHb3w9zLj8Czoluop451HlYa5/r2OfHtA+TPhzNbGFwoFUOYvLnqauMMenAU8D3rLU1rvN0lzHmAqDcWrvUGHO64zhHKhGYCFxnrV1ojLkTmA381G2s7jPGZNM+ajQMqAKeMMZcaa39V6Ta9FwxttYetLgaYx6i/doLwBPE+LDPYc7lO8DTHcV3kTEmRPs8qRXRytcdBzsXY8wE2n+gV3asxFMALDPGTLHWlh3oe1w61P8TAGPMV4ELgDNj9YPRIXRp2lqvMMb4aS/ED1trn3adp4emARd1LEMbADKNMf+y1l7pOFdPFAFF1trOEYonaS/GXnQW8Im1tgLAGPM0cDIQsWIcb8PUJcBpHV+fAXzsMMuRepb2c8AYMwpIwoMTr1trV1lr+1trC621hbT/wk6MxUJ8OMaYc4GbgIustQ2u8/RAV6a29QTT/snub8A6a+3+y1Z5hLX2ZmttQcfvxuXAPI8WYjp+p7cZY47ueOtM9l5q10u2AlONMakdP2tnEuGb0TzXMz6MbwF3dlxwb+LTFaK86AHgAWPMaqAF+KoHe2Lx5s9AMvBqRy9/gbX2areRuu5gU9s6jtVT04AvA6uMMSs63vtRx2yB4s51wMMdH/Y24dGpkTuG2Z8EltF+SWo5EZ6NSzNwiYiIOBZvw9QiIiKeo2IsIiLimIqxiIiIYyrGIiIijqkYi4iIOKZiLCIi4piKsYiIiGMqxiIiIo79fygHINIgAi6AAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = 50\n",
    "x_min = -8.0\n",
    "x_max = 8.0\n",
    "\n",
    "fig = plt.figure(figsize=(8, 6))\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
}