Skip to content

Instantly share code, notes, and snippets.

@jcboyd
Created January 15, 2023 15:55
Show Gist options
  • Save jcboyd/78e1ba35ed5af22a5c641b3f1b96e149 to your computer and use it in GitHub Desktop.
Save jcboyd/78e1ba35ed5af22a5c641b3f1b96e149 to your computer and use it in GitHub Desktop.
Metropolis-Hastings algorithm
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"id": "e3b6c103",
"metadata": {},
"source": [
"# Metropolis-Hastings\n",
"\n",
"The acceptance probability of sample $x' \\sim g(x' | x)$ is given by,\n",
"\n",
"$$A(x, x') = \\min(1, \\frac{f(x')}{f(x)}\\cdot\\frac{g(x | x')}{g(x' | x)}),$$\n",
"\n",
"where $f(x)\\propto p(x)$, and $g(x | x')$ governs the transitions."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "d7dd7ba7",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f269a33f430>"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA5wAAAEvCAYAAAApc84jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhtUlEQVR4nO3deVhUZR/G8e8AgqgImgouuJSWkiWuRHuJmWVlmllZLq9puZRJbmjuC6amprm3aVnappWVZZZWRppbi1tZGW7glqioIDDvH0+AJCLgDGeGuT/XdS7OzJyZuVEYzu88m81ut9sRERERERERcTAvqwOIiIiIiIhI8aSCU0RERERERJxCBaeIiIiIiIg4hQpOERERERERcQofqwOIiIiIiIjkJiMjg9TUVKtjyAWUKFECb2/vPI9RwSkiIiIiIi4nNTWVv/76i4yMDKujSB6CgoIICQnBZrPl+rgKThERERERcSl2u50DBw7g7e1NaGgoXl4aCehq7HY7p06d4uDBgwBUrlw51+NUcIqIiIiIiEtJS0vj1KlTVKlShVKlSlkdRy7A398fgIMHD1KpUqVcu9fqUoGIiIiIiLiU9PR0AHx9fS1OIheTeUHg7NmzuT6uglNERERERFzShcYFiuu42P+RCk4RERERERFxChWcIiIiIiIixViXLl1o06aNJe+tSYNERERERMQtFHUPW7u9aN+vOFILp4iIiIiIiDiFCk4REREREREHee+997jmmmvw9/fnsssuIyoqiuTkZH788UdatGhBhQoVCAwM5JZbbmHTpk05nmuz2Zg7dy6tW7emVKlS1KtXj7i4OHbt2sWtt95K6dKluf766/njjz+ynjNy5EjCw8OZO3cuoaGhlCpVigcffJCkpKQLZszIyCA2NpZatWrh7+9PgwYNeO+997Ie/+eff+jYsSMVK1bE39+fOnXq8NprrxXq30MFp4iIiIiIiAMcOHCAhx9+mP/9739s376d1atX07ZtW+x2OydOnKBz58589913/PDDD9SpU4e77rqLEydO5HiNMWPG0KlTJ7Zs2ULdunV55JFHeOKJJ4iJiWHDhg3Y7Xb69OmT4zm7du3inXfe4eOPP2bFihVs3ryZXr16XTBnbGwsCxcuZM6cOWzdupV+/frx6KOPsmbNGgCGDRvGtm3b+Oyzz9i+fTuzZ8+mQoUKhfo30RhOERERERERBzhw4ABpaWm0bduWGjVqAHDNNdcAcPvtt+c4dt68eQQFBbFmzRpat26ddX/Xrl158MEHARg0aBCRkZEMGzaMli1bAtC3b1+6du2a47XOnDnDwoULqVq1KgAzZszg7rvv5oUXXiAkJCTHsSkpKYwfP54vv/ySyMhIAC6//HK+++475s6dyy233EJ8fDwNGzakSZMmANSsWbPQ/yZq4RQREREREXGABg0a0Lx5c6655hrat2/P/Pnz+eeffwBITEyke/fu1KlTh8DAQMqWLcvJkyeJj4/P8RrXXntt1n5wcDCQXbRm3nfmzBmOHz+edV/16tWzik2AyMhIMjIy2Llz53kZd+3axalTp2jRogVlypTJ2hYuXJjVVbdnz54sXryY8PBwBg4cyPfff1/ofxO1cIqIiIiIiDiAt7c3K1eu5Pvvv+eLL75gxowZDB06lHXr1tGzZ0+OHDnCiy++SI0aNfDz8yMyMpLU1NQcr1GiRImsfdu/0/Lmdl9GRkahMp48eRKATz75JEeRCuDn5wdAq1at+Pvvv/n0009ZuXIlzZs3p3fv3kyePLnA76cWThEREREREQex2WzccMMNjBo1is2bN+Pr68vSpUtZu3YtTz/9NHfddRdXX301fn5+HD582CHvGR8fz/79+7Nu//DDD3h5eXHVVVedd2xYWBh+fn7Ex8dTu3btHFtoaGjWcRUrVqRz5868+eabTJs2jXnz5hUqm1o4RUREREREHGDdunWsWrWKO+64g0qVKrFu3ToOHTpEvXr1qFOnDm+88QZNmjTh+PHjDBgwAH9/f4e8b8mSJencuTOTJ0/m+PHjPP300zz44IPnjd8ECAgIoH///vTr14+MjAxuvPFGkpKSWLt2LWXLlqVz584MHz6cxo0bc/XVV5OSksLy5cupV69eobKp4BQREREREXGAsmXL8s033zBt2jSOHz9OjRo1eOGFF2jVqhUhISH06NGDRo0aERoayvjx4+nfv79D3rd27dq0bduWu+66i6NHj9K6dWtmzZp1wePHjBlDxYoViY2N5c8//yQoKIhGjRoxZMgQAHx9fYmJiWH37t34+/tz0003sXjx4kJls9ntdnuhnikiIiIiIuIEZ86c4a+//qJWrVqULFnS6jgubeTIkSxbtowtW7ZY8v4X+7/SGE4RERERERFxChWcIiIiIiIi4hTqUisiIiIiIi5FXWrdh7rUioiIiIiIiCVUcIqIiIiIiIhTqOAUERERERERp1DBKSIiIiIiIk6hglNEREREREScQgWniIiIiIiIOIUKThEREREREYvUrFmTadOmZd1OSEigRYsWlC5dmqCgIMtyOYqP1QFERERERETEmDp1KgcOHGDLli0EBgZaHeeSqeAUERERERFxEX/88QeNGzemTp06VkdxCHWpFRERERERcYBbb72VPn360KdPHwIDA6lQoQLDhg3DbrcDcPDgQe655x78/f2pVasWixYtyvH8mjVr8v7777Nw4UJsNhtdunSx4LtwLLds4UxLS2Pz5s0EBwfj5aWaWURERESkOElLSyMtLY3U1FRzvm+3w6lTuR5rt9tJS0ujVKlS2Gw2x4cpVQoK8LoLFiygW7durF+/ng0bNtCjRw+qV69O9+7d6dKlC/v37+frr7+mRIkSPP300xw8eDDruT/++COdOnWibNmyvPjii/j7+zv++ylibllwbt68mWbNmlkdQ0REREREnKBGjRrMmTOHlJQUALxOn6bRzTdf8Hg/Z4Y5eRJKl8734aGhoUydOhWbzcZVV13FL7/8wtSpU7nlllv47LPPWL9+PU2bNgXglVdeoV69elnPrVixIn5+fvj7+xMSEuLwb8UKbllwBgcHA7B+/XoqV65scRoREREREXGktLQ0Tpw4QY0aNShZsiQkJ1sdKd+uu+66HC2tkZGRvPDCC2zfvh0fHx8aN26c9VjdunWLxUy0eXHLgjOzG23lypWpVq2axWlERERERMSRzpw5w+nTp/H19cXX1xdKlDAtjblITU3l119/pX79+uZYRytVyvGv6UHcsuAUEREREREPYrNduFtriRJk+Pubx51RcBbQunXrctz+4YcfqFOnDnXr1iUtLY2NGzdmdanduXMnx44dsyBl0dGMOyIiIiIiIg4SHx9PdHQ0O3fu5O2332bGjBn07duXq666ijvvvJMnnniCdevWsXHjRh5//PFiMTFQXlRwioiIiIiIOEinTp04ffo0zZo1o3fv3vTt25cePXoA8Nprr1GlShVuueUW2rZtS48ePahUqZLFiZ1LXWpFREREREQcpESJEkybNo3Zs2ef91hISAjLly/Pcd9jjz2W4/ayZcucGa/IqeAUERERERGPc/DgQRISEjh79iylSpWievXqlL7AONHDhw+ze/fuHPfZbLYcM85K7lRwioiIiIiIRzl69Ch79uyhRo0alC5dmsTERH777Tfq169PiRIlcn2Ot7c39evXL+Kk7k8Fp4i4p/37YfduqFMHKla0Oo2IWOHECdi6FSpUgCuuMLNYiojkQ2JiIhUrVqRChQoA1KhRg6SkJA4fPkzlypUv+LwLFaOZVq9e7ciYxYIKThFxffv3w4YNsHFj9paQkP141arQsGHOrUYNnXyKFCeJibB5c85t167sxwMDoVEjaNIEGjc2m4pQEY+Snp5Oenp61m2bzYaX1/lzpGZkZJCcnExISEiOY8uWLUtycnKer//zzz8DUKpUKapWrVrsZ5h1BBWcIlLk8jr/s9v/3Vm7FiZNgnXrchaXmby8oHJl2Lcvezt3EH5QEFx3HQwbBtdf78j4IlJU/vgDRoyAr76CAwdyPyYkBP75B5KS4OuvzZYpMNAUnk88Ae3bq/gUcUP2rBODi9u6dWuO21WqVKFKlSrnHZeWlgac31pZokQJzpw5k+trlyxZkpo1a1KqVCnS09NJSEhgx44dXH311fi6wNqfVrrY/5EKThFxLfv2wcCB8NZb2fd5eUFYWHarRePG0KCBWeD5xAn46SfT2rFli/n6669w7BisWGG2hx+GCROgenWrvisRKYikJBg3DqZNg7NnzX02G1x1FYSHZ/dkCA83XerPnjVda8/tBfHTT+Z1vvrKbDNnwvTp5rNDRFyet7c3AKmpqfluRfxv8Wdz4EWmMmXKUKZMmazbpUuXZuvWrRw6dIiqVas67H3c0alTp4ALdzdWwSkiLsGPM0QzBa4aD8nJ5uSyWzf43//MCWKpUrk/MSAAbrzRbJlSU83J58yZ8Oqr8PbbsGwZDBhgitkLzEAnIhZLTze/s889BwcPmvtatoQhQ8yFpgv87tp8SwDh/27dAPDhLGFs46cRS2HiRPjmG9PltkcPGDPGjPsUEZfl4+NDqVKlOHToECVKlMi1a2ym1NRUAM6ePZuvIjMjIwOA5ORkfHyyy6EzZ87g5eV1wVbO//Lz8+PUqVP5Pr64sdvtnDp1ioMHDxIUFJR1keC/bPaCtFO7iL179xIaGsqePXuoVq2a1XFEpIBy/i2wcx8f8gLPcgV/mruuv960RDhiqvHNm+GZZ8zJJkCVKqa1s2NH03IqIq7h66/N7+q/46O46iqYMgVatbpoV9iLdtOPjzcXm5YsMXeWKwejRkHPnuCja+8irio1NZW//vorq0C8kLS0NPbt20fVqlVzFJB5OXDgAH5+fpQvXx4wxdO+ffsICAggMDDwos+32+0cOHCAkiVLZr2GpwoKCiIkJOSCxb4KThEpcpmfR/XYxov0pQVfArCPKlRdNMl0gXXkWCu7HT74APr3NzPbAjRrBi++aMZ5ioh1/vzT/G4uXWpuBwXByJHQqxdcZDbITPkaFw6wZg307Wu62wLUr28+B26/vTDJRaQIZGRkZLVgXkhCQgK33norq1evzjERUF4+/fRTBg8ezKhRo7j22mtZsGABK1as4NNPP6VChQoMGjSI4OBgoqOjAZg5cyYNGjSgRo0aHD9+nFdffZUvv/yS999/n9q1a1/y9+muSpQoccGWzUy6rCcilujAYhbSCV/OkoIvk+lPLDGcfKTMxZ9cUDYbtGsHd99txoSNGwfr15tuuK+8Ap07O/49ReTi1q41v5dJSeDtDU8+aVoeL7vMOe93yy1mfOf8+TB0qBnv3by56cI7erQmFRJxQV5eXpQsWTLPY3x8fPj777/x8fG56LGZ2rZty/79+xk8eDAJCQmEh4czf/78rMasn376iZo1a2a93t69e5k+fToJCQmUK1eOxo0bs3jxYq3LmQ9q4RSRItfDNo85PIkXdj7hLp5iBn9xOfCf1ghnOXDAtHK8+665PXWq6conIkVnxQpo2xZOn4aICHPx5+qrC/VS+W7hPNfRo2YW61mzzO2nnjIXpNTVXsTtqDZwbfpUFZGiNXky83gCL+zM5knu4eOsYrPIVK4MixdDv37mdr9+5sTT/a6/ibinxYvhnntMsdmqlZlFtpDFZqGVL28mFps1y1SsM2aYScr+XS5BREQcQwWniBQNu910WxswAIAJDKIXs7Bb9THk5QUvvABjx5rbY8dC795wkYkJROQSzZ4NjzxiCruHHzYzSF9oFuqi0LMnLFxouvQuWAAdOkBKinV5RESKGRWcIuJ8GRmmC+u4cQAMJpYYJgAWj5ey2cw4rswWjtmzzey1F5mcQEQKwW43nwG9epn9Xr3gzTfBFRZMf/RReP99k+WDD+Dee83yTCIicslUcIqIc6WlmW5qM2aY2zNn8jyDrc30Xz17wltvmeURFi+GNm3g30WMRcQBMjLg2WdNLwcwXdhfesm1xkvedx98+qlZ6/OLL8z6n8eOWZ1KRMTtudAnvYgUOykppnvaggWmu9obb5hWDVf00EPw8cfg7w+ffQZ33KGTTRFHyLzoNHWquT11quvOCNu8OaxcaZZmWbsWbrsNDh60OpWIiFtTwSkiznH2rGkx+OAD003tvfdMtzVXduedOU82mzdXS6fIpbDbzbJDmRedFixw/RmhIyNh9WqoVAm2bIGbb1bRKSJyCVRwiohzREfD55+byUA++cR0U3UHN9xgFoevWBE2bTItM5q9VqRwJk7M7q7+/vvQqZPVifKnQQP49lsIDYWdO+HBB81FNBERKTAVnCLieAsXmvFZAG+/DVFR1uYpqGuvNSfHPj6wZIk5aRaRgvnsM4iJMfvTp5seD+7kyivNRbMyZcxFqIEDrU4kIuKWVHCKiGNt2gRPPGH2hw83sz26o5tuMifJYE6aV6ywNo+IO/n9d7Pkid0O3bvDk09anahw6tUzF9AApk0zrbUiIlIgKjhFxHEOH4b774czZ+Duu2HECKsTXZonnzQny3a7mVTo99+tTiTi+o4fN62ZSUlw/fWmt4MrThCUX/ffD0OGmP3HHzfjOkVEJN8KVHDGxsbStGlTAgICqFSpEm3atGHnzp05jjlz5gy9e/fmsssuo0yZMrRr147ExMQcx8THx3P33XdTqlQpKlWqxIABA0hLS7v070ZErJOWZmakjY+H2rXN+nqutORBYdhsZjmXyEhz8nzffXDihNWpRFxXRoYZp7l9O1Spkr22pbsbPdosk3L6tClAjxyxOpGIiNso0NngmjVr6N27Nz/88AMrV67k7Nmz3HHHHSSfszhyv379+Pjjj3n33XdZs2YN+/fvp23btlmPp6enc/fdd5Oamsr333/PggULeP311xk+fLjjvisRKXoxMfDVV2YNu2XLsJULwmYj182t+PmZk+YqVcxJ9GOPmZNqETnf6NHw4Yfm92bpUggJsTqRY3h7m+60l18Ou3fDI49AerrVqURE3ILNbi/89IuHDh2iUqVKrFmzhptvvpmkpCQqVqzIW2+9xQMPPADAjh07qFevHnFxcVx33XV89tlntG7dmv379xMcHAzAnDlzGDRoEIcOHcI3H1dC9+7dS2hoKHv27KFatWqFjS8ijrJkielyCvDOO9C+faELy0uZEDav97zkiWbXrzfjOlNTTVfhkSMv8QVFipmlSyHzAvPrr5vlUIqIU3/3z/Xzz6bHw6lTMHgwxMY68MVFpLBUG7i2S+rvlpSUBED58uUB2LhxI2fPniXqnBkp69atS/Xq1YmLiwMgLi6Oa665JqvYBGjZsiXHjx9n69atub5PSkoKx48fz9pOqEubiOv45RezdAjAoEHQvr3T3upCLaZF0nLarBnMnWv2R42CZcuc/IYibmTr1uwlT/r2LdJis0hdey288orZnzDBrC8sIiJ5KnTBmZGRwTPPPMMNN9xA/fr1AUhISMDX15egoKAcxwYHB5OQkJB1zLnFZubjmY/lJjY2lsDAwKwtLCyssLFFxJH++cesr3nqFLRoAePGWZ3Iubp0gaefNvuPPWZOskU8XebnwMmTcPvtMHmy1Ymc66GHzDrDYD4T9DkgIpKnQhecvXv35tdff2Xx4sWOzJOrmJgYkpKSsrZt27Y5/T1F5CLsdtOi8eefULOmWW/T29vqVM43eTLcdps5uW7f3szIK+KpMpc92bXLfA4sWWLWry3unn/efA4kJ5tJhE6etDqRiIjLKlTB2adPH5YvX87XX3+do590SEgIqampHDt2LMfxiYmJhPw7cUBISMh5s9Zm3g65wOQCfn5+lC1bNmsLCAgoTGwRcaQFC2D5cjM5yAcfwGWXWZ2oaJQoYU6qg4PNJEIayyme7J13zKRaPj7ma4UKVicqGj4+5nOgWjWzXNLQoVYnEhFxWQUqOO12O3369GHp0qV89dVX1KpVK8fjjRs3pkSJEqxatSrrvp07dxIfH09kZCQAkZGR/PLLLxw8eDDrmJUrV1K2bFl1lRVxFwcOQL9+Zn/0aGjY0No8Ra1ixezxnJMmwbp11uYRsUJiIvTubfaHDoVGjazNU9QqVoSXXzb7M2bA2rXW5hERcVEFKjh79+7Nm2++yVtvvUVAQAAJCQkkJCRw+vRpAAIDA+nWrRvR0dF8/fXXbNy4ka5duxIZGcl1110HwB133EFYWBiPPfYYP/30E59//jnPPfccvXv3xs/Pz/HfoYg4lt0OPXvCsWPQpEn2WCZPc9990LGjWSKlSxd1rRXPYrdDr15mPcoGDWDIkHw/1bKJv5yhZUvo2tX8e/zvf2adThERyaFABefs2bNJSkri1ltvpXLlylnbkiVLso6ZOnUqrVu3pl27dtx8882EhITwwQcfZD3u7e3N8uXL8fb2JjIykkcffZROnToxevRox31XIuI877xj1tkrUQJefdUzxmtdyPTpZp3BHTvMUikinmLJEtOV3sfHLIGSjyXNiq0XXoDKleG338wM1iIiksMlrcNpFa21I2KRQ4cgLAwOHzZjF/MospyxDueltIA47ZPuo49Ma6eXl+lS929vDpFiKzERrr7atG5e5HMgN85YM7PI1uG8kHM/B9atM70/RKTIqDZwbZe0DqeIeJi+fU2xec01EBNjdRrXcO+98Oijpmtt167qWivFW2aX+iNHIDy8QF1pi7V77zXLpWRkmK61qalWJxIRcRkqOEUkfz78MHvpk9dec1oXOrcc3/Xii9lda4cPtzqNiPMsXgxLl2Z3pS1RwupErmP6dDNL7y+/QGys1WlERFyGCk4Rubhjx0yrBkD//tC4saVxCsOphWz58jBvntl/4QX44QcHvKiIi0lIgD59zP6wYWayIMlWsSK89JLZHzvWFJ4iIqKCU0Ty4dlnzVIoV12VY7yWW7ZGOss998Bjj2XPWqvZKqU4yexKe/SoWQZJXepz9+CD0KYNpKWZrrVpaVYnEhGxnApOEcnbF1+Y2WhtNnjlFfD3tzqR63rxRTNb5c6d6lorxcvbb8OyZaYLrbrSXpjNBrNmQVAQbNgAU6ZYnUhExHIqOEXkwk6cgO7dzf7TT8MNN1ibx9WVK5eza21cnLV5RBwhIQGeesrsDxsG115rbR5XV7kyTJ1q9ocPNxegREQ8mApOEbmwoUMhPp4/qUXpF8ep22x+tG4NnTqZLohPPqkudeL+BgzI7ko7eLDVadxD587QsiWkpMDjjxfR2iwiIq5JBaeI5O6XX2DmTAB6MI9TlLY4kPM4fCzqlClmIqGff4Y5cxyaVaRIffcdvPmm+WWYP19dafPLZjO9HUqXNv+Gb79tdSIREcuo4BSR89nt8MwzkJHBuzzAKqKsTuReLrvMzFIJpgvioUPW5hEpjPT07K60jz/ulrNTW6p69ex1SgcNguRka/OIiFhEBaeInO/DD+Grr8DPjwFMsjqNe+rRA8LDzZIyQ4danUak4ObNgy1bzNjk8eOtTuOeoqOhZk3YuxcmTrQ6jYiIJVRwikhOKSlmGRSA/v35m5qWxnFb3t7Za/K9/LKZsVLEXRw+nH2hZMwYqFDB2jzuqmRJmDzZ7E+cCH//bW0eEclh5syZ1KxZk5IlSxIREcH69evz9bzFixdjs9lo06aNcwMWEyo4RSSnadPgzz+hShVNEHKpbrjBrM1pt0OfPmaNThF38Nxz8M8/0KABPPGE1WncW9u2cOutcOYMDBxodRoR+deSJUuIjo5mxIgRbNq0iQYNGtCyZUsOHjyY5/N2795N//79uemmm4ooqftTwSki2Q4cyB57+PzzUKaMtXmKg+efh4AAWLcOFiywOo3IxW3cmL28z4wZ4ONjbZ5/OXxyr6J6P5vNXMjz8oJ33oFvv3VGXBEpoClTptC9e3e6du1KWFgYc+bMoVSpUrz66qsXfE56ejodO3Zk1KhRXH755UWY1r2p4BSRbEOGwMmTEBEBjzxidZrioXJlsxYfmBbjpCRr84jkJSPDTBRkt5vPgCK+gl/URWWRadAge03jvn3NhEwi4nAnTpzg+PHjWVtKSkqux6WmprJx40aiorInRfTy8iIqKoq4PNbQHj16NJUqVaJbt24Oz16cqeAUEePHH+H1183+iy+aq/HiGE8/DXXrwsGDMHKk1WlELuzNNyEuzvRumKQJwxxqzBgIDITNm+G116xOI1IshYWFERgYmLXFxsbmetzhw4dJT08nODg4x/3BwcEkJCTk+pzvvvuOV155hfnz5zs8d3GnM0oRMa0Zffua/U6dTAunOI6vL0yfbvZnzIBff7U2j0hukpKyxxgOG2bGcYvjVKyYfcFp6FD1dhBxgm3btpGUlJS1xcTEOOR1T5w4wWOPPcb8+fOpoEnUCkwFp4iYRcnj4swi5Re4GuipHNbFr0ULuP9+05Xu6adNkS/iSkaPhsREuPJKsw6vOF7v3nDVVaa3Q+Z4eRFxmICAAMqWLZu1+fn55XpchQoV8Pb2JjExMcf9iYmJhISEnHf8H3/8we7du7nnnnvw8fHBx8eHhQsX8tFHH+Hj48Mff/zhlO+nuFDBKeLpkpOzWzWGDFGrhjNNmWKWSfj6a3jvPavTiGTbti27FX76dNMqL45XogRMnWr2X3wRfv/d2jwiHsrX15fGjRuzatWqrPsyMjJYtWoVkZGR5x1ft25dfvnlF7Zs2ZK13Xvvvdx2221s2bKF0NDQoozvdlRwini6iRNh3z6zOHl0tNVpio1cW0Vr1WTEmX+XmomONsW+iNXsdtPqnpYG990HLVtanah4a9UK7roLzp7NXvNYRIpcdHQ08+fPZ8GCBWzfvp2ePXuSnJxM165dAejUqVNWl9ySJUtSv379HFtQUBABAQHUr18fX12ky5MKThFP9vffpuAEszh5yZLW5vEAExloivu9e7P/7UWs9NFHsGoV+PmZVnhxvilTzHIzH38Mn39udRoRj9ShQwcmT57M8OHDCQ8PZ8uWLaxYsSJrIqH4+HgOHDhgccriwWa3u99Aor179xIaGsqePXuoVq2a1XFE3FfHjvDWW2ZR8q++ynVgotsvR+BEeX165vXvZn/3PWjfHkqVMl3q1I1ZrHL2LNSvD7/9ZrrUjxvn9Lcs6s+Uwp7l5Pk77Igzp+ho07326qvhp5/A29sBLyrimVQbuDa1cIp4qi1bTLEJ5mq7Ksui064dREbCqVMwYoTVacSTzZ9vis2KFWHQIKvTeJbhw6FcOdi6FRYtsjqNiIjTqOAU8VRDhpivDz0EDRtam8XT2GymCzPAq6+aE06Ronb8ePYyHSNGQNmylsaxgsNmoS6MoCAY/O+Y7uHD4QIL1IuIuDsVnCKeaM0a+OwzM4ZozBir03im66+Htm0hIyN7lmCRojRxIhw6ZJZB6dHDoS9taSHnTp56ynSp//tvmDPH6jQiIk6hglPE09jt2VfVu3eH2rWtzePJJkwwRf+nn5oxtCJFZd++7AmCJkwwy3VI0fP3z25lHjvWtDqLiBQzKjhFPM1HH8EPP5gJa4YNszqNZ6tTB5580uz3729aO0WKwrBhcPo03HADtGljdRqncvnW1q5dTSvz4cOaJVhEiiUVnCKeJD09e+zmM89A5cqWxhHM2K2yZWHz5uxJnESc6eef4fXXzf7kyS5UeXkoH5/s2YFfeAEOHrQ2j4iIg6ngFPEkb7wB27ZB+fIaN+gADmk5qVgxu4vz0KFw5oxTsopkGTjQdK1v3x6uu87qNAJm5uomTeDkySJZmkZEpCip4BTxFGfOZC/BERMDgYHW5pFszzwD1apBfDxMn251GinOVq6Ezz83YzZjY61OI5lsNjOWFmD2bNi929I4IiKOpIJTxFPMnm0KmmrVoHfvHA+5/Bin4s7f30wYAjB+PBw5Ym0eKZ7S02HAALPfqxdccYW1eSSn5s0hKgrOnjVd7UVEigkVnCKe4Pjx7G5aI0eaAkdcy6OPQoMGkJSkpWrEOd58E376yfRu0IRhrimz1fnNN+GXX6zNIiLiICo4RTzB5Mmm1axuXejc2eo0khtvb5g0yezPmgV//GFtHileTp+G554z+0OGwGWXAerd4HKaNDFja+327AneRETcnApOkeIuMTF7qv1x48yMiOKaWrSAli1Nl7qYGKvTSHEybRrs3QvVq8PTT1udRvIydqy5ALV8OXz3ndVpREQumQpOkeJu7FhIToZmzeD++61OIxczcaJpXnr3XbNeqsilOnw4e0KaceOgZElr80jerrwSunUz+4MHm9ZOERE3poJTpDj780+YO9fsT5igfnLu4NproUsXs6+TTXGEcePMOO6GDeGRR6xOI/kxfLi5MLB2rWnpFBFxYyo4RYqz0aNN98w77oDbbrM6jeTXqFHg5wdr1sCKFVanEXf2999mTDCYi05e+rPvFqpWhb59zf6wYZCRYW0eEZFLoL88IsXVjh3wxhtmP3PJDXEPoaHQp4/Zj4nRyaYU3vDhkJrKKm7H1rKFJgZyJwMGQECAmVn4/fetTiMiUmgqOEWKq1GjTKFy773QtKnVaeQc+ZoZNCYGypY1J5tvv21ZVnFjv/ySddFpMBMAVZhu5bLLoF8/sz9ihFlHVUTEDangFCmOfv0Vliwx+6NHW5tFCueyy2DQILM/bBikplqbR9zPkCFgt/MuD7ABXXQqKg5daqZfPwgKgu3bYfFiR0cVESkSBS44v/nmG+655x6qVKmCzWZj2bJlOR7v0qULNpstx3bnnXfmOObo0aN07NiRsmXLEhQURLdu3Th58uQlfSMico4RI8xkM+3bQ4MGVqeRwurbF0JC4K+/sid/EsmP774zk814ezOUcVankcIKCjJda8H0WklLszSOiEhhFLjgTE5OpkGDBsycOfOCx9x5550cOHAga3v7P93BOnbsyNatW1m5ciXLly/nm2++oUePHgVPLyLn27wZPvjAXE4fOdLqNHIpSpc2Fw8AxoyBEyeszSPuwW7Pbh3v1o3fudLaPHJpnn4aKlSA33/PHpcvIuJGClxwtmrVirFjx3J/Huv5+fn5ERISkrWVK1cu67Ht27ezYsUKXn75ZSIiIrjxxhuZMWMGixcvZv/+/YX7LkQk2/Dh5usjj0BYmLVZ5NJ16wZ16sChQzBlitVpxB18/DF8/z34+2dfsCgEh3YNlXzJ9d87oAz9D/97AWH0aHWvFxG345QxnKtXr6ZSpUpcddVV9OzZkyNHjmQ9FhcXR1BQEE2aNMm6LyoqCi8vL9atW5fr66WkpHD8+PGs7YSu8ovkbt26rG50l3KiKdY572TTtwQP/v7vLMOTJ8PBg9YGFNeWnm4mnALTJbtKFWvziEPMopfpXr97N7z2mtVxREQKxOEF55133snChQtZtWoVzz//PGvWrKFVq1ak/zu7WkJCApUqVcrxHB8fH8qXL09CQkKurxkbG0tgYGDWFqZWG5HcDRtmvnbqZFrFpFh4jwfYQGM4eRLGaTye5GHhQti2DcqVy+5WK27vNKWyLySMHQtnzlgbSESkABxecD700EPce++9XHPNNbRp04bly5fz448/snr16kK/ZkxMDElJSVnbtm3bHBdYpLj49ltYuRJ8fLILTykW7Hj9u6wFMHu2mURI5L/OnMnu2RATYyackeKjRw+oVg327oX5861OIyKSb05fFuXyyy+nQoUK7Nq1C4CQkBAO/qdLWFpaGkePHiUkJCTX1/Dz86Ns2bJZW0BAgLNji7gXuz27yOzWDWrVOu8Qjcdyb6uIgqgoOHs2e5yuyLlmzoQ9e0xR0qeP1WnE0UqWhKFDzf748XDqlLV5RETyyekF5969ezly5AiVK1cGIDIykmPHjrFx48asY7766isyMjKIiIhwdhyR4umrr2DNGlLwJXTuUBWVxdWEf1s5Fy2Cn3+2Nou4lqQkU4SAWT7D39/aPOIc//sf1KwJCQmmt4OIiBsocMF58uRJtmzZwpYtWwD466+/2LJlC/Hx8Zw8eZIBAwbwww8/sHv3blatWsV9991H7dq1admyJQD16tXjzjvvpHv37qxfv561a9fSp08fHnroIapocgORgjundXMuT7CXUIsDidM0bgwdOpj/88zxXCIAEyfC0aNQr54Zwy3Fk69vdm+WCRPMuG4RERdX4IJzw4YNNGzYkIYNGwIQHR1Nw4YNGT58ON7e3vz888/ce++9XHnllXTr1o3GjRvz7bff4ufnl/UaixYtom7dujRv3py77rqLG2+8kXnz5jnuuxLxJCtWQFwc+PsTi4qQYm/sWDNO99NPYc0aq9OIKzhwAKZONfvjx5ufDyfLq4u+elQ4WadOULs2HD4MM2ZYnUZE5KJsdrvdbnWIgtq7dy+hoaHs2bOHatWqWR1HxDp2OzRtChs3Qv/+2CZPsjqROFHWp3Xv3jBrFkREmIsNOsP3bE88AfPmQWQkrF173s+DfjwcL68zJ2f9e+d4z0WL4NFHzWzEf/0FgYHOeVMRN6HawLU5fQyniDjRRx+ZYrN0aRg40Oo0UlSGDTP/5+vWwQcfWJ1GrLRjB7zyitl//nlVl57ioYcgLAz++Se7dVtExEWp4BRxVxkZ2UsgPP00VKxobR4pOiEh8OyzZj8mxsxcK55p6FBIT4d77oGbbrI6jRQVb28YOdLsT51qCk8RERelglPEXS1dCj/9BAEB0L+/1WmkqPXvby4y/P57dguXeJa4ONPC7eUFsbFWp/EoLjF+tV07uOYaOH4cpkwpwjcWKT5mzpxJzZo1KVmyJBEREaxfv/6Cx37wwQc0adKEoKAgSpcuTXh4OG+88UYRpnVfKjhF3NG5rZv9+kH58tbmkaIXEJC9HufIkZqt0tPY7TBokNnv0gWuvtrSOGIBLy+zBA7AtGlw5IilcUTczZIlS4iOjmbEiBFs2rSJBg0a0LJlSw4ePJjr8eXLl2fo0KHExcXx888/07VrV7p27crnn39exMndjyYNEnEzNhs8yBKW8BDHCKQmu0kiyOpYUgTO+7ROTTXjuP74A0aPzl4uQYq/5ctNN9qSJU0rdx5/CzWss3jI9WzNbjfLJW3eDIMHq6VbPFZhaoOIiAiaNm3KSy+9BEBGRgahoaE89dRTDB48OF+v0ahRI+6++27GjBlT6OyeQC2cIm7Gi3RGMhKAF3hWxaYn8/WFcePM/sSJcOiQtXmkaKSnm+ICoG/fPItNKeZstuxWzhkz4AItMyKe4sSJExw/fjxrS0lJyfW41NRUNm7cSFRUVNZ9Xl5eREVFERcXd9H3sdvtrFq1ip07d3LzzTc7LH9xpYJTxM08xGLqsYOjlONF+lodR6zWvr1p4Th5EnSF1TMsXAhbt5olMTK71Yrnat0amjSB5GSYpKWxxLOFhYURGBiYtcVeoNX/8OHDpKenExwcnOP+4OBgEhISLvj6SUlJlClTBl9fX+6++25mzJhBixYtHPo9FEcqOEXcSVoaIzBXsycxgBOUtTiQWM7LyyyHATBnjuleK8XX6dPZY3eHDDFFp3g2m810qQeYORPyOFkWKe62bdtGUlJS1hYTE+PQ1w8ICGDLli38+OOPjBs3jujoaFavXu3Q9yiOVHCKuJNFi7iS3zlEBV6ij9VpxFU0bw4tW5rlUTSOs3ibMQP27oXQUOijzwD51513QmSkuSAxYYLVaUQsExAQQNmyZbM2Pz+/XI+rUKEC3t7eJCYm5rg/MTGRkJCQC76+l5cXtWvXJjw8nGeffZYHHnjggq2okk0Fp4i7OHs26yr2RAZykgCLA4lLmTDBtHS8/TZs3Gh1GnGGo0ezJ4UZM8ZMGCQCOVs558yBffuszSPi4nx9fWncuDGrVq3Kui8jI4NVq1YRGRmZ79fJyMi44DhRyaaCU8RdLFwIf/5JIpWYRS+r04irCQ+Hjh3Nfj5n1xM3M2ECHDtm1l589FGr04irad4cbroJUlI0W61IPkRHRzN//nwWLFjA9u3b6dmzJ8nJyXTt2hWATp065eiSGxsby8qVK/nzzz/Zvn07L7zwAm+88QaP6vP4onysDiAi+ZCamjUhzAQGc4rSFgcSlzRmDLzzDnz5JXzxBdxxh9WJxFHi42H6dLM/YQJ4e+d4WEufSFYr5223wfz5MHAgVK9udSoRl9WhQwcOHTrE8OHDSUhIIDw8nBUrVmRNJBQfH4+XV3bbXHJyMr169WLv3r34+/tTt25d3nzzTTp06GDVt+A2tA6niDuYOxeefBJCQvBP+JMz+FudSCyQr0/rfv3MIvDh4bBhw3mFibipLl1gwQK45Rb4+uvzKkwVnMVfvs/Wbr/d/Iw88YTpXiviAVQbuDZ1qRVxdSkpMHas2R8yRMWm5G3oUAgMhC1bTDdscX+bNmX/X06cqOpS8pa5Lucrr8Du3ZZGEREBFZwiru/ll82slFWrQvfuVqcRV1ehAjz3nNkfMsSszynuy243rdZ2OzzyCDRrZnUicXU33QQtWkBaWvbFShERC6ngFHFlp0/D+PFmf+hQzUop+fPUU3DFFWY9vsw1OsU9LV0K33xjfvc1EYzkV2Yr5+uvw65dlkYREVHBKeLKZs2C/fvNxA//+5/VacRd+PmZrpcAkyebCWfE/aSkwIABZr9/f00AI/kXGQl33QXp6dnFp4iIRVRwiriqEyeyF/AeMcIUESL5df/9ZoKZM2fgnGndxY3MmAF//gmVK8OgQVanEXfz78zmLFoEW7dam0VEPJoKThEL2Gx5bwC8+CIcPgxXXgmdOlmaV9yQzQZTppivb70F69ZZnUgK4tCh7IJh3DgoU8baPOJ+GjWCdu3M+N/hw61OIyIeTAWniCv65x/TFRJMdygfLZkrhdCoEXTubPYzJ54R9zBiBBw/Dg0bZv8fihTUqFHmotMHH8DGjVanEREPpYJTxBVNmgRJSXDNNfDgg1anEXc2bhyULg1xcfDOO1ankfzYutWsvQswdSr8u/D4RXtFiPzX1VdDx45mP3P2ahGRIqaCU8TFVCLRdKcF06XOS7+mcgmqVMke/zdokBnTKa6tf3/IyMgehytyKUaONL1kVqyA776zOo2IeCCdyYq4mBhi4dQps97evfdaHUeKg2efhWrV4O+/TYuZuK4VK8xWokT2TMMil+KKK7JnOX/uOXWtF5Eip4JTxIVUYw89mW1ujB2rvnLiGKVKZc94PH68WZ9TXE9amrk4APD001C7trV5pPh47jnw9YU1a+DLL61OIyIeRgWniAt5jrH4kWq60UVFWR1HipOHHzat5idPwrBhVqeR3MybB9u2QYUKGm8njhUaCj17mv2hQ9XKKSJFSgWniIu4nD/4H6+aG2rdFEfz8sruTvvKK/DTT9bmkZyOHcteumLUKAgKsjKNFEcxMaa3w48/wscfW51GRDyICk4RFzGSkZQgjc+4E2680eo4Uhxdf72Z9dhuh7591crhSkaNgiNHoF496NHD6jRSHAUHm997ML0cMjKszSMiHkMFp4gLCGMrHVkEmG61Ik4zcSL4+5uxXG++aXUaAdi8GaZPN/vTpmndXXGeAQMgMBB+/lnLJIlIkVHBKeICRjMcL+y8T1s20djqOFKc1aiR3XXz2Wfh6FFr83i6jAwzti4jw7Q+33GH1YnERTlkHdZy5bInphoxwkxUJSLiZCo4RSzWiI204wMysDGc0VbHEU8QHQ1hYXDoEAwZYnUazzZ/PqxbBwEBWrJGisYzz5iJqX77DRYutDqNiHgAFZwiFhuDmTF0ER3ZxtUWpxGP4OsLs/9dfmfuXIiLszaPp0pMhMGDzf64cVClirV5xDMEBGT/3I0aBSkp1uYRkWJPBaeIhW7la+7iM87iw0hGWh1HPMnNN0OXLma/Z091rbNC//5mdtpGjaBXL6vTiBsrcHfbXr3MBY74eJg1q0iziojnUcEpYhEbGUxkIABzeYI/uSL7MUeM1RG5mIkToXx5s0TKjBlWp/EsX39tJm2y2WDOHPD2tjqReBJ/f9O6CWYZrmPHLI0jIsWbCk4Ri7TnXZqygROUYTTDrY4jnqhiRVN0glkmYc8ea/N4CD9bCjtu7wnAS/Ze2Jo11QUlKXpduphleI4eheeftzqNiBRjKjhFLFCCVMZjJmuZxAAOUcniROKxunY163MmJ5vJRMTpBjCJuuwkgWAtgyTW8fHJLjSnTYO9ey2NIyLFlwpOEQs8wVyu4E8SCGYK0VbHkWKg0N2wvbyyu3R+8AF88kmR5PVYf/zBUMYB0I+pJBFkbR7xbK1bw003wZkzZpkUEREnUMEpUtSOH89a/mQEo0imjMWBxONdcw3062f2+/SBU6eszVNc2e3Qpw/+nOFLmrOYh6xOJJ7OZsvuVv/66/Drr5bGEZHiSQWnSFGbOJGKHGYHV/EK3axOI2KMGAGhobB7t5lERBzvvfdgxQpS8KUXswAN2hQXcN110K4dZGRkL5ciIuJAKjhFitL+/TBlCgAxxJKOj8WBRP5Vpkz2TLWTJ8O2bdbmKW6OH88aIxtLDL9zpbV5RM41frzpVv/JJ7BmjdVpRKSYKXDB+c0333DPPfdQpUoVbDYby5Yty/G43W5n+PDhVK5cGX9/f6Kiovj9999zHHP06FE6duxI2bJlCQoKolu3bpw8efKSvhERtzByJJw+zVquZxltrE4jktN998E998DZs2YyobNnrU5UfERHmwtOtWszAbUiiYu58kp44gmzP3Cg6f4tIuIgBS44k5OTadCgATNnzsz18YkTJzJ9+nTmzJnDunXrKF26NC1btuTMmTNZx3Ts2JGtW7eycuVKli9fzjfffEOPHj0K/12IuIPt2+GVVwAYyETUnU5c0qxZEBQE69fDuHFWpykeli0zv/s2G7z8MimUtDqRyPmGDzc9HdavN92/RUQcxX4JAPvSpUuzbmdkZNhDQkLskyZNyrrv2LFjdj8/P/vbb79tt9vt9m3bttkB+48//ph1zGeffWa32Wz2ffv25et99+zZYwfse/bsuZT4IkXrvvvsdrDb77vPbi4fa9NWsC0vhX1ert5+2zzR29tuj4srzE+7x8v8tw/mgP0gFex2sE9goNP+j7Vpu9BWICNHmifVrm23p6YW+udfpKipNnBtDh3D+ddff5GQkEBUVFTWfYGBgURERBAXFwdAXFwcQUFBNGnSJOuYqKgovLy8WLduXa6vm5KSwvHjx7O2EydOODK2iPN99x18+KFZgiI21uo0Inl76CF45BFIT4dHHwUNeSgkO6/yPypymC00yJqdWsRlPfssBAfDrl0wb57VaUSkmHBowZmQkABAcHBwjvuDg4OzHktISKBSpZyL3Pv4+FC+fPmsY/4rNjaWwMDArC0sLMyRsUWcy26HAQPMfrduUK+etXlE8mPmTDNr7R9/mPGHUmBPMoe7+Iwz+NGRRaTiZ3UkkbyVKZO9HueoUWayKxGRS+QWs9TGxMSQlJSUtW3T7IniTpYuhR9+gFKlzKRBIkXMZst7y1VQECxYYA6YPx8++qgoI7u9K9nJCzwLwGAmsI2rLU4kkk+PP24mETp0iNGBk/P/mSHihmbOnEnNmjUpWbIkERERrF+//oLHzp8/n5tuuoly5cpRrlw5oqKi8jxesjm04AwJCQEgMTExx/2JiYlZj4WEhHDw4MEcj6elpXH06NGsY/7Lz8+PsmXLZm0BAQGOjC3iPCkpMGiQ2Y+OhipVrM0jbq3AReOluu227NbNxx+H/3y2ywWcPcubPEopTvMlzZnO01YnEsm/EiWyhn70ZzLV2GNxIBHnWLJkCdHR0YwYMYJNmzbRoEEDWrZseV6dkmn16tU8/PDDfP3118TFxREaGsodd9zBvn37iji5+3FowVmrVi1CQkJYtWpV1n3Hjx9n3bp1REZGAhAZGcmxY8fYuHFj1jFfffUVGRkZREREODKOiPWmTTNjYUJCzFTzIu5m3Di45ho4dMgUnXa71Ylc35gxNGUD/xBEF17H7h6diUSy3X8/a7iZUpxmIvrbJcXTlClT6N69O127diUsLIw5c+ZQqlQpXn311VyPX7RoEb169SI8PJy6devy8ssvk5GRkaPukdwV+K/gyZMn2bJlC1u2bAHMREFbtmwhPj4em83GM888w9ixY/noo4/45Zdf6NSpE1WqVKFNmzYA1KtXjzvvvJPu3buzfv161q5dS58+fXjooYeootYfKU7274cxY8z+88+DWubFHfn5waJF4OsLy5eb7rVyYXFxWcvJPMkc9lHN4kAihWCz0ZcXSceLh1nMTXxjdSKRfDlx4kSOiUZTUlJyPS41NZWNGzfmmOjUy8uLqKiorIlOL+bUqVOcPXuW8uXLOyR7cVbggnPDhg00bNiQhg0bAhAdHU3Dhg0ZPnw4AAMHDuSpp56iR48eNG3alJMnT7JixQpKlsxed2zRokXUrVuX5s2bc9ddd3HjjTcyT7OhSXEzeDAkJ0NEhJnpU8RdXXMNjB9v9vv1g99/tzaPqzp5Eh57DDIyeJOOvEMHqxOJFNpPhDOf7gBM52m8SLc4kcjFhYWF5ZhoNPYCKwMcPnyY9PT0PCc6vZhBgwZRpUqVHEWr5M5mt7tf/6i9e/cSGhrKnj17qFZNV4/FBcXFwfXXm/3166Fp0xwPa9IFcSX5+iuQkQFRUfD11+YiynffgY+P07O5le7d4eWXITSUoD0/k0RQgV8ir/8LfW5IYRT2LM9mg8s4zO/UoRzHeJLZzOXJS3pNEWfJrA22bdtG1apVs+738/PDz+/8GcL3799P1apV+f7777OG/YFpOFuzZs0Fl2rMNGHCBCZOnMjq1au59tprHfeNFFMaWCLiaBkZ8PS/k4R07XpesSnilry8zKy1gYGwbl32Uj9iLFxoik2bDRYuLFSxKeJqjlAha/3YsTxHEP9YnEgkbwEBATkmGs2t2ASoUKEC3t7eeU50eiGTJ09mwoQJfPHFFyo280kFp4ijvf46bNhgxmxmdkMUKQ5CQ+GVV8z+tGkaz5lp7VrTugkwZAjcemuhX6rIZyIWuYjZ9ORXrqYCRxjFCKvjiDiEr68vjRs3zjHhT+YEQOe2eP7XxIkTGTNmDCtWrKBJkyZFEbVYUMEp4khJSRATY/ZHjDCz04oUJ+3amQXhAXr1Ml1sPdnu3XD//ZCaCm3bwujRVicScah0fOjLiwD0YhZX86vFiUQcIzo6mvnz57NgwQK2b99Oz549SU5OpmvXrgB06tSJmMxzOuD5559n2LBhvPrqq9SsWZOEhAQSEhI4efKkVd+C21DBKeJIo0fDwYNw1VXw1FNWpxFxjmHD4KGHIC3NFKCeOonQ8eNwzz1myZhGjUy3Wi/9WZXi5yua8wH340M6L9JXgzilWOjQoQOTJ09m+PDhhIeHs2XLFlasWJE1kVB8fDwHDhzIOn727NmkpqbywAMPULly5axt8uTJVn0LbkOTBok4yo4dZjbPtDT49FNo1eqCh6p7nLiSQv0VOH0abrvNjOe86iozUVa5cg7P5rLS0+G+++CTT6ByZTM52Dl/j/Q7Lq7iUiYNOldN/mI79ShJCrz3nrnYJOIiVBu4Nl2KFXEEux2eecYUm61b51lsihQL/v6wbJkZ17lzJzz4IJw9a3WqojNwoCk2S5aEDz/MUWyKFEe7qcUk/p0s7NlnzUUnEZF8UMEp4gjLl8Pnn0OJEjBlitVpRAqk0BPVhITARx9B6dLw5ZdmjU5P8PLL2b/nCxZoJmrxGBMYzB6qwd9/g7oRikg+qeAUuVQpKdkn2v36QZ061uYRKUrh4fDmm6Y6nTnTbMXZ6tXQs6fZHzXKtOyKeIhTlGYAk8yN2FiIj7c2kIi4BRWcIpfqhRfgjz9Ma89zz1mdRqTotWljTj4B+vaFL76wNI7T7Nplxq2lpcHDD5vJk0Q8zBI6wE03mS61zz5rdRwRcQMqOEUuxe+/Zy+DMGmSWXtTxBMNHAidO5vJdB580EwmVJzs3Qt33w1Hj0JEhFmPVDMDiUeywfTp4O1tJg/66COrA4mIi1PBKVJYGRlmsfeUFGjRAjp2tDqRiHVsNpg7F26+2axH27w5rFxpdSrH2LkTbrgBfvsNqlc3kyX5+1udSsQ64eEQHW32e/UySwSJiFyACk6Rwnr1VVizBkqVMifaau0QT+fnZ2ZubdECkpNNi+B771md6tJs2mS6D8bHm+Vfvv3WdJ8X8XQjR8Lll8O+fRATY3UaEXFhKjhFCuPAAejf3+yPGQO1almbR8RVlCkDH38M7dubZVIefBDmzbM6VeGsXg233gqHDkHjxqbYrF7d6lQirqFUqezf7VmzYO1aa/OIiMtSwSlSGE8/bboNNmli9kUkm58fvP02PPGEWaP2iSfMpEKFXYHeCsuWwZ13wokTcNtt8NVXULGi1alEXEvz5tC1q9nPHGIiIvIfKjhFCmrZMtNN0NvbrMfn45PrYYVe21CkOPD2htmzYehQc3vIENMrICPD2lz58dprZjbalBQzA++nn0LZslanEnE4h/ydmjwZKlWC7duzZ6sWETmHCk6RgkhKgt69zf6AAdCggbV5RFyZzQZjx8KUKeb2lCnwv/+ZZUVc1QsvmIwZGebru+9CyZJWpxJxXeXLw4wZZn/8eNi61do8IuJyVHCKFMTgwbB/P9SuDcOHW51GxD306wevv25aPRcsgPvug4QEq1PldOqUWUM0c2x2//559mAQcRdF0tumfXu45x4zbrt7d/foySAiRUYFp0h+ffcdzJlj9ufN07IIIgXRuTN88IEZ3/npp1Cvnvk9coUT088/h/r1zdqCAM8/b9bVVf93kfyx2WDmTDNpWFyc6U4vIvIvFZwi+XHmjLlqC9Ctm5lERMQDOLR15N57zclo48Zw7JiZTOjmm2HbNkfHzp/ERHjkETM50F9/QWiomWF34MB8PV3jtEXOERoKEyaY/cGDYc8ea/OIiMtQwSmSH+PHw44dEBxsWj5EpHAaNoQffoCpU6F0abOUQng4DBtmLuwUhYwM0122bl0zm66Xl+n2u20btG6d41AVlSIF0LMnREbCyZPQq5d7zUwtIk6jglPkYn7+OXvmvZdegnLlrM0j4u58fOCZZ8yslvfea8Z9jR0L115rlh9xpu3bzdqa3bubVtZGjWD9ejOhUZkyzn1vkeLOy8tczClRApYvh8WLrU4kIi5ABadIXk6fNl3u0tLMRCft2lmdSKT4CA01ywy9/z5UqQK//27W9Wvd2kwudPSoY97n9Gn48EMzjrRBA/j2W9O6OmUKrFtnuviKiGOEhWUvh9SzJ8THW5tHRCynglMkLwMHwtatJBBMpQ/nYfOyqUudiCPZbNC2renO2ru3uf3JJ9Cli1nbLyoKZs0ys0MXxLFjsGgRPPAAVKhg1tNcuNC0prZubZZu6NdPs9CKOMOQIRARYZYSe/RRSE+3OpGIWMhmt7tfB/u9e/cSGhrKnj17qFatmtVxpLhavtxM8w7cyWd8zp3nHZLXb48KUinunPLXY9s2eOcdM6PtL7/kfOy660zhWL36hZ9/9KiZ+Oerr0xxmal6dbj/flOA3nBDvn9B9XssUjBZnwt//GHGZ588CWPGwHPPWRlLijnVBq5NBadIbg4cMOPJDh9mCv14lim5HqaCUzyZ0/967NoFS5eaLS6u4M+vV8+0nt5/vxmrWYhfSv0eixRMjs+FhQtNV3Zvb7O02HXXWZZLijfVBq5NBafIf2VkQKtW8MUX0KABfj+tIxW/XA9VwSmerEj/euzfb8ZhrlgByckXPq5ECbPUyv33m1loL5F+j0UKJsfngt0OHTua2aBr1YItW6BsWauiSTGm2sC1afCKyH+9+KIpNv394a23SL0692JTRIpQlSpmApKePa1OIiL5ZbPB7Nmmh8Jff0GfPqbVU0Q8iiYNEjnX5s0waJDZnzrVzLYnIiIihRMYCG++aZZMeeMNM5mXiHgUFZwimU6dMkugnD1rlkDp0cPqRCIe59xZoP+7iYibuuEGGD7c7PfsaVo7RcRjqOAUyRQdDTt2mK57L7+sM1wRERFHGTrUFJ4nTphxnWlpVicSkSKiglMEzCyYc+eaInPhQrNun4iIiDiGj4/pWlu2rBnTOWaM1YlEpIio4BT5+294/HGzP2AANG9ubR4REZHiqGZNc3EXYOxYWL3ayjQiUkRUcIpnS0424zWPHoUmTXTFVaQIaJymiAd76CGzNmdGBjzwgMZzingAFZziuTIyzB+9n36C4GD44APw9bU6lYiISPE2ezY0bgxHjpiLvidPWp1IRJxIBad4rrFj4f33zULxH3wAoaFWJxJxK+7SUukuOUU8hr8/LFtmLvb+8gt06mQuAotIsaSCUzzTBx/AiBFmf84cuP56a/OIiIh4kmrVzIR9vr7m66hRVicSESdRwSme56ef4LHHzP4zz8D//mdpHBEREY8UGQnz5pn90aPh3XetzSMeZ+bMmdSsWZOSJUsSERHB+vXrL3js1q1badeuHTVr1sRmszFt2rSiC+rmVHCKZzl0yIwXOXUKWrSASZOsTiQiIuK5Onc262Bn7m/ebG0e8RhLliwhOjqaESNGsGnTJho0aEDLli05ePBgrsefOnWKyy+/nAkTJhASElLEad2bCk7xHKmpZka8v/+G2rVhyRKzLpiIOJzGTYp4pkL97j//PLRsCadPm4vCiYlFllc815QpU+jevTtdu3YlLCyMOXPmUKpUKV599dVcj2/atCmTJk3ioYcews/Pr4jTujeHF5wjR47EZrPl2OrWrZv1+JkzZ+jduzeXXXYZZcqUoV27diTqg0WczW6Hp56Cb74xi05/9BGUK2d1KhEREfHxgcWL4corYc8eaNfOXCQWKaATJ05w/PjxrC0lJSXX41JTU9m4cSNRUVFZ93l5eREVFUVcXFxRxfUYTmnhvPrqqzlw4EDW9t1332U91q9fPz7++GPeffdd1qxZw/79+2nbtq0zYohkmz3bjBOx2eDtt6FePasTiYiISKagIHMxODAQ1q6FXr3MxWKRAggLCyMwMDBri42NzfW4w4cPk56eTnBwcI77g4ODSUhIKIqoHsUp/Ql9fHxy7duclJTEK6+8wltvvcXtt98OwGuvvUa9evX44YcfuO6665wRRzzde++Z1k2ACRPgrruszSMiLkNdfEVcyFVXmYvCd98Nr7wCVatq9lopkG3btlG1atWs2+r66hqc0sL5+++/U6VKFS6//HI6duxIfHw8ABs3buTs2bM5mq/r1q1L9erV82y+TklJydE8fuLECWfEluLok0/g4YfN+l7dusGAAQ59eY1TExERcaBWrWD6dLM/ejRMnGhtHnErAQEBlC1bNmu7UMFZoUIFvL29zxvWl5iYqAmBnMDhBWdERASvv/46K1asYPbs2fz111/cdNNNnDhxgoSEBHx9fQkKCsrxnIs1X8fGxuZoHg8LC3N0bCmOvvrKjANJS4OHHoK5cy9YCapwFHFv+h0WKUb69IHx483+oEEwc6a1eaTY8fX1pXHjxqxatSrrvoyMDFatWkVkZKSFyYonh3epbdWqVdb+tddeS0REBDVq1OCdd97B39+/UK8ZExNDdOaU2cC+fftUdErevv8e7r0XUlLMjHcLF4K3t9WpREREJD9iYuDkSVN49ukDpUpB165Wp5JiJDo6ms6dO9OkSROaNWvGtGnTSE5Opuu/P2edOnWiatWqWeNAU1NT2bZtW9b+vn372LJlC2XKlKF27dqWfR/uwOlrQgQFBXHllVeya9cuWrRoQWpqKseOHcvRynmx5ms/P78cTeLHjx93ZmRxd5s2mS45yclwxx1m+ZMSJaxOJSIiIgUxdqz5W/7ii/D446bo7NDB6lRSTHTo0IFDhw4xfPhwEhISCA8PZ8WKFVkTCcXHx+Plld0ZdP/+/TRs2DDr9uTJk5k8eTK33HILq1evLur4bsXpBefJkyf5448/eOyxx2jcuDElSpRg1apVtGvXDoCdO3cSHx+v5mtxjF9/NUXm8eNw002wdClowLiIiIj7sdlg6lRTdL78Mjz6qCk677nH6mRSTPTp04c+ffrk+th/i8iaNWti18zJheLwMZz9+/dnzZo17N69m++//577778fb29vHn74YQIDA+nWrRvR0dF8/fXXbNy4ka5duxIZGakZauXS/f47REXBkSPQrBksX27+MImIiIh7stlgzhx45BEzJ8MDD8DKlVanEpECcHgL5969e3n44Yc5cuQIFStW5MYbb+SHH36gYsWKAEydOhUvLy/atWtHSkoKLVu2ZNasWY6OIZ7m77+heXNITIRrr4XPPoOyZXMcoslDRNybfodFiq+8fr/tdm9YsABOnzY9l+67Dz7/3PRkEhGXZ7O7Ydvw3r17CQ0NZc+ePVSrVs3qOGK17dvNmM2//zZreH3zDVSqdN5hOlkVERGxTl5nnHkXnP/upKRAmzawYgWUKWPW2W7Z0pERxU2pNnBtTlmHU6TIfPstXH+9KTbr1IFVq3ItNkVERMTN+fnB+++bHk0nT8Ldd8Orr1qdSkQuQgWnuK933jFjNo8dg8hIsxRK1apWpxIRERFnKVUKPv0UOnaE9HTo1g1Gjsy7+VRELKWCU9yP3Q5Tppip0VNT4f77TctmhQpWJxMREZELsNkuvBWIry+88YZZqxNg1ChTeJ496/DMInLpVHCKe0lPh2eegWefNbefegrefRf8/S2NJSIiIkXIZoPx480Mtl5e8NprZrmUEyesTiYi/6GCU9zH6dPw4IMwfbq5PXmyWQza29vaXCIiImKNJ56ADz80XW0//xxuuQUOHLA6lYicQwWnuIcjR8x4zQ8+MF1pFi82rZyaelZERMSztW4Nq1dDxYqwebOZ12HbNqtTici/VHCK6/v2W2jY0EwKFBRkFnzu0MHqVCIiIuIqmjaFuDgzY/3ff0NEhBnnKSKWU8EpristzUwEcOutsGcPXHklrF0LN998wac4bEICERERcS9XXGEuTt96q1k2pVMneOwxjesUsZgKTnFNe/bA7bebqc4zMqBLF9i4EcLCrE4mIiIirqpCBfjySxgzxszx8OabppfUhg1WJxPxWCo4xfUsXQoNGpiutAEBsGiRmX2uTBmrk4mIiIir8/aG556DNWugenX44w8zrnPyZHMRW0SKlApOcR2nT0OvXtC2LfzzDzRrBlu2wCOPWJ1MRERE3M0NN5jziAceMMN0BgyAVq0gIcHqZCIeRQWnuIYNG8yA/9mzze1Bg+C77+Dyy887VOM0RUREJF/KlYN33oF588ya3V98YXpRLVsGdrvV6UQ8ggpOsdahQ9C9u2nN3LoVgoPNH4MJE6BECavTiYiIiLuz2cy5xoYNcM01cPAg3H+/ae3cudPqdCLFngpOsUZaGrz0kpl59uWXzVXGRx+Fn3+GFi2sTiciIiLFTVgYrFsHMTFmTe/PPzcF6KBBmslWxIlUcErR++YbaNwYnnoKjh2D8HDTffaNN6BSJavTiYiISHHl7w/jx8Ovv8Ldd8PZszBxItStC2+9pW62Ik6gglOKzr59ZgKgW24xLZnly5sxmxs2mIH9IiIiIkWhTh1Yvhw+/tis37l/P3TsaM5RfvrJ6nQixYoKTnG+ffvg2Wfhqqvg7bfNWIonn4TffjNfvb2tTigiIiKeqHVr09o5bhyUKmWWZGvUyAzz+fVXq9OJFAsqOMV5fvsNHn8catWCKVMgORmuv960aM6eDZddZnVCERER8XQlS8KQIbBjBzz4oFmrc9EiM77znnvg+++tTiji1lRwiuNt2gTt25vxEK+8YsZH3HwzfPaZGavZqJHVCUVERERyCg2FJUtg40ZzHmOzmW63N9yQfR6jMZ4iBaaCUxwjIwNWrYKWLc2EQO+9Zz6U77kH1q6FNWvgzju1WKaIiIi4tkaNzNqdO3aYnlolSpiutnfdBQ0bmuFBKSlWpxRxGyo45dL8/jsMGwaXXw5RUWYNTW/v7CVOPvrIdKMtIJvtwpuIiIiI0115JcyfD3/9ZeaiKF3aTCj0yCNQpQr07g3r16vVU+QiVHBKwf3zD8ydawrJK6+EsWPh77+hbFnz4fv772aJk2uusTqpiIiIuLm8LkI760J0jveoVhXbC5MpnxzPcEaZYvPoUZg1CyIizPqesbGwZ49zwoi4OZvd7n6XZfbu3UtoaCh79uyhWrVqVsfxDCdOwJdfmm4kH32U3ZXEywvuuAM6d4b77jPrWzmAWjJFREQkP5xxJpvXeYg9Ld0MI1qwAJYuhdOns590++1meZW77oLgYMcHk1ypNnBtPlYHEBdlt5vpwD/7LHuyn7S07Mfr1zdFZseOULmydTlFREREipK3t7nYfscdcPy4mbdi4UIzX8WqVWYDMxa0VSuzRUSAj067xTOphVOyHT0Kq1ebAnPFCti7N+fjtWub9ao6dYLwcKc2Q6qFU0RERC5VYc9y82zhvNBr/vUXvPkmLFtmZuw/V7ly0KKFKT6jokDnrw6l2sC1qeD0VKmpZlKfdevM9sMPZuzluUqWhNtuy746V7u2QyOoqBQRERFnckbBma/3S0iAzz83F/G/+MLMf3GuqlVNq2fm1qSJmZRICkW1gWtTwekJkpJg+3az/fKLKTA3bsx9Su+6dc3SJq1amTWnHDQmMzcqOEVERMSZLCs4z5WWBj/+mD1MafNmSE/PeYy3txmuFBFhepGFhUG9elCxok6Y8kG1gWtTZ/LiIjUV9u2D3buzi8vMbf/+3J9Trpz5YLvuOvO1WTMoX75Ab3uxz0D3u5whIiIi4kA+PhAZabbRoyE52Vz4z+xltm6dGcb0009mO1f58qbwPHe7/HIIDYVSpaz5fkQKSAWnq7PbzQyxBw/CoUOQmGg+lOLjzVIk8fFmO3Ag7+quShXzIRUWBk2bmiKzdm1dNRMREREpSqVLm15kN9+cfd++fabwXL/eTNq4fbsZE3r0KKxda7b/qlABatSA6tWzt9BQMztupUpmCwoyKwpIrmbOnMmkSZNISEigQYMGzJgxg2bNml3w+HfffZdhw4axe/du6tSpw/PPP89dd91VhIndkwpOZ7PbzXTZyclmO3XKFJDHjpmurrl9PXLEFJiZRWZuXV9zU7Kk+bC58srsq2BhYaabbGCg877HPKieFREREbmIqlWhbVuzZTp9Gn77DbZty+61tmOH6c128iQcPmy2jRsv/Lre3qZbbqVK5mvFiqYIDQzM/WvZsqbltHRps5UqVWwL1iVLlhAdHc2cOXOIiIhg2rRptGzZkp07d1KpUqXzjv/+++95+OGHiY2NpXXr1rz11lu0adOGTZs2Ub9+fQu+A/fh1mM497/5JpUvu8wUdRkZ5389d0tPP/92Wpr5mtv+2bN5bykpZjtz5vyvZ87kLDAdoUyZ7A+MatVyXs3K3Czo56+CUkRERNxRXmfADh3D6Wh2u2mkyOzldu62Z092o8WxY455v5IlswtQf3/w8zP35fa1RIm8Nx8fUwRnfv3vvre3KXDP3c69z2YzW+b+v18PHDlClY4dCzSGMyIigqZNm/LSSy8BkJGRQWhoKE899RSDBw8+7/gOHTqQnJzM8uXLs+677rrrCA8PZ86cOY75ty6m3LqFs/Sjj1odoWAyf2HLlLnw1aXAQNNFIvNKVOZX9dMXERERcRi3vWhus5nzxqAguPbaCx+Xmmp6yh06lLPnXFLShXvZHT9uGkvObTDJbEw5csSp39alyJzf98SJExw/fjzrfj8/P/z8/M47PjU1lY0bNxITE5N1n5eXF1FRUcTFxeX6HnFxcURHR+e4r2XLlixbtuyS8xd3bl1wpoeFmasp/73aYbPlfoUk8yqJzXb+FZb/XmnJvApzoSs0F7qyk7mfeSUosztCqVLmtUVEREREnM3X13TVrVq14M/NyDBdek+dyu61l5ycXXz+t6df5n5qat49BP/bs/C/X/PqoZieblp3M7dzejWmp6bC9u2EhYXl+DZGjBjByJEjz/v2Dh8+THp6OsHBwTnuDw4OZseOHbn+kyQkJOR6fEJCQsH/fT2MWxecyZ9/TjlNfex0bnsFUEREREQKzssru+GkYkWr01xU8t69EBrKtm3bqHpOgZ1b66YUPbcuOEVERERERAACAgIoW7bsRY+rUKEC3t7eJCYm5rg/MTGRkJCQXJ8TEhJSoOMlW/GcdkpERERExIVkjvrKbZOi5evrS+PGjVm1alXWfRkZGaxatYrIyMhcnxMZGZnjeICVK1de8HjJphZOERERERHxKNHR0XTu3JkmTZrQrFkzpk2bRnJyMl27dgWgU6dOVK1aldjYWAD69u3LLbfcwgsvvMDdd9/N4sWL2bBhA/PmzbPy23ALKjjdUF5Xwgo7HbeuromIiIhcmsKeTznj3E7y1qFDBw4dOsTw4cNJSEggPDycFStWZE0MFB8fj9c5a5Bef/31vPXWWzz33HMMGTKEOnXqsGzZMq3BmQ9uvQ5nQdbacTdFvf6TCk4RERER9+J+Z/HO4Qm1gTuzdAznzJkzqVmzJiVLliQiIoL169dbGadY0PgAERERERFxFZYVnEuWLCE6OpoRI0awadMmGjRoQMuWLTl48KBVkYo9FaMiIiIiniGv8z6d+0lRsqxLbUREBE2bNuWll14CzMxQoaGhPPXUUwwePDjHsSkpKaSkpGTd3rNnD/Xr12f9+vVUrlw53+8ZGnrhx/bsKVh+R7xuXs8TEREREXGWSzn3dTUHDhygWbNm/P3331SvXt3qOPIflkwalJqaysaNG4mJicm6z8vLi6ioKOLi4s47PjY2llGjRp13f7NmzRyWyVnFn4pKEREREXE1xfEcdc+ePSo4XZAlBefhw4dJT0/PmgUqU3BwMDt27Djv+JiYGKKjo7NuHz16lFq1avHrr78SGBjo9Lziek6cOEFYWBjbtm0jICDA6jhiAf0MCOjnQPQzIPoZEEhKSqJ+/frUq1fP6iiSC7dYFsXPzw8/P7/z7g8NDaVs2bIWJBKrHT9+HICqVavqZ8BD6WdAQD8Hop8B0c+AkPX/7uPjFqWNx7Fk0qAKFSrg7e1NYmJijvsTExMJCQmxIpKIiIiIiIg4mCUFp6+vL40bN2bVqlVZ92VkZLBq1SoiIyOtiCQiIiIiIiIOZlm7c3R0NJ07d6ZJkyY0a9aMadOmkZycTNeuXS/6XD8/P0aMGJFrN1vxDPoZEP0MCOjnQPQzIPoZEP0MuDrLlkUBeOmll5g0aRIJCQmEh4czffp0IiIirIojIiIiIiIiDmRpwSkiIiIiIiLFlyVjOEVERERERKT4U8EpIiIiIiIiTqGCU0RERERERJxCBaeIiIiIiIg4RbEoOD/55BMiIiLw9/enXLlytGnTxupIYpGUlBTCw8Ox2Wxs2bLF6jhSRHbv3k23bt2oVasW/v7+XHHFFYwYMYLU1FSro4kTzZw5k5o1a1KyZEkiIiJYv3691ZGkiMTGxtK0aVMCAgKoVKkSbdq0YefOnVbHEgtNmDABm83GM888Y3UUKWL79u3j0Ucf5bLLLsPf359rrrmGDRs2WB1LzuH2Bef777/PY489RteuXfnpp59Yu3YtjzzyiNWxxCIDBw6kSpUqVseQIrZjxw4yMjKYO3cuW7duZerUqcyZM4chQ4ZYHU2cZMmSJURHRzNixAg2bdpEgwYNaNmyJQcPHrQ6mhSBNWvW0Lt3b3744QdWrlzJ2bNnueOOO0hOTrY6mljgxx9/ZO7cuVx77bVWR5Ei9s8//3DDDTdQokQJPvvsM7Zt28YLL7xAuXLlrI4m53DrZVHS0tKoWbMmo0aNolu3blbHEYt99tlnREdH8/7773P11VezefNmwsPDrY4lFpk0aRKzZ8/mzz//tDqKOEFERARNmzblpZdeAiAjI4PQ0FCeeuopBg8ebHE6KWqHDh2iUqVKrFmzhptvvtnqOFKETp48SaNGjZg1axZjx44lPDycadOmWR1LisjgwYNZu3Yt3377rdVRJA9u3cK5adMm9u3bh5eXFw0bNqRy5cq0atWKX3/91epoUsQSExPp3r07b7zxBqVKlbI6jriApKQkypcvb3UMcYLU1FQ2btxIVFRU1n1eXl5ERUURFxdnYTKxSlJSEoB+5z1Q7969ufvuu3N8Hojn+Oijj2jSpAnt27enUqVKNGzYkPnz51sdS/7DrQvOzJaLkSNH8txzz7F8+XLKlSvHrbfeytGjRy1OJ0XFbrfTpUsXnnzySZo0aWJ1HHEBu3btYsaMGTzxxBNWRxEnOHz4MOnp6QQHB+e4Pzg4mISEBItSiVUyMjJ45plnuOGGG6hfv77VcaQILV68mE2bNhEbG2t1FLHIn3/+yezZs6lTpw6ff/45PXv25Omnn2bBggVWR5NzuGTBOXjwYGw2W55b5pgtgKFDh9KuXTsaN27Ma6+9hs1m491337X4u5BLld+fgxkzZnDixAliYmKsjiwOlt+fgXPt27ePO++8k/bt29O9e3eLkotIUenduze//vorixcvtjqKFKE9e/bQt29fFi1aRMmSJa2OIxbJyMigUaNGjB8/noYNG9KjRw+6d+/OnDlzrI4m5/CxOkBunn32Wbp06ZLnMZdffjkHDhwAICwsLOt+Pz8/Lr/8cuLj450ZUYpAfn8OvvrqK+Li4vDz88vxWJMmTejYsaOucrmx/P4MZNq/fz+33XYb119/PfPmzXNyOrFKhQoV8Pb2JjExMcf9iYmJhISEWJRKrNCnTx+WL1/ON998Q7Vq1ayOI0Vo48aNHDx4kEaNGmXdl56ezjfffMNLL71ESkoK3t7eFiaUolC5cuUcdQBAvXr1eP/99y1KJLlxyYKzYsWKVKxY8aLHNW7cGD8/P3bu3MmNN94IwNmzZ9m9ezc1atRwdkxxsvz+HEyfPp2xY8dm3d6/fz8tW7ZkyZIlREREODOiOFl+fwbAtGzedtttWT0dvLxcsgOHOICvry+NGzdm1apVWctgZWRksGrVKvr06WNtOCkSdrudp556iqVLl7J69Wpq1apldSQpYs2bN+eXX37JcV/Xrl2pW7cugwYNUrHpIW644YbzlkT67bffVAe4GJcsOPOrbNmyPPnkk4wYMYLQ0FBq1KjBpEmTAGjfvr3F6aSoVK9ePcftMmXKAHDFFVfoireH2LdvH7feeis1atRg8uTJHDp0KOsxtXgVT9HR0XTu3JkmTZrQrFkzpk2bRnJyMl27drU6mhSB3r1789Zbb/Hhhx8SEBCQNXY3MDAQf39/i9NJUQgICDhvzG7p0qW57LLLNJbXg/Tr14/rr7+e8ePH8+CDD7J+/XrmzZunXk4uxq0LTjBLH/j4+PDYY49x+vRpIiIi+Oqrr7T+jogHWblyJbt27WLXrl3nXWRw45WfJA8dOnTg0KFDDB8+nISEBMLDw1mxYsV5EwlJ8TR79mwAbr311hz3v/baaxfthi8ixUfTpk1ZunQpMTExjB49mlq1ajFt2jQ6duxodTQ5h1uvwykiIiIiIiKuS4OcRERERERExClUcIqIiIiIiIhTqOAUERERERERp1DBKSIiIiIiIk6hglNEREREREScQgWniIiIiIiIOIUKThEREREREXEKFZwiIiIiIiLiFCo4RURERERExClUcIqIiIiIiIhTqOAUERERERERp/g/yMxjAwUlgT4AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1000x300 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def norm(mu, sigma):\n",
" return lambda x : np.exp(-0.5 * (x - mu) ** 2 / sigma)\n",
"\n",
"f = lambda x : 0.5 * (norm(-2, 1)(x) + norm(+2, 1)(x))\n",
"g = lambda x, y : norm(y, 1)(x)\n",
"\n",
"iters = 10000\n",
"samples = []\n",
"\n",
"x0 = np.random.randn()\n",
"\n",
"for t in range(iters):\n",
" x = x0 + np.random.randn() \n",
" A = min(1, f(x) / f(x0) * g(x, x0) / g(x0, x))\n",
"\n",
" if np.random.rand() <= A:\n",
" x0 = x\n",
"\n",
" samples.append(x0)\n",
"\n",
"samples = samples[iters // 10:] # apply burn-in\n",
"\n",
"xmin, xmax = -6, 6.\n",
"xs = np.linspace(xmin, xmax, num=100)\n",
"\n",
"fig, ax = plt.subplots(figsize=(10, 3))\n",
"_ = ax.hist(samples, bins=100, color='blue', label='samples')\n",
"ax.set_xlim([xmin, xmax])\n",
"\n",
"ax2 = ax.twinx()\n",
"ax2.plot(xs, f(xs), color='red', label='pdf')\n",
"\n",
"fig.legend()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment