{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Manual definition of Feynman-Kac models\n", "\n", "It is not particularly difficult to define manually your own `FeynmanKac` classes. Consider the following problem: we would like to approximate the probability that $X_t \\in [a,b]$ for all $0\\leq t < T$, where $(X_t)$ is a random walk: $X_0\\sim N(0,1)$, and \n", "$$ X_t | X_{t-1} = x_{t-1} \\sim N(x_{t-1}, 1).$$\n", "\n", "This probability, at time $t$, equals $L_t$, the normalising constant of the following Feynman-Kac sequence of distributions: \n", "\n", "\\begin{equation}\n", "\\mathbb{Q}_t(dx_{0:t}) = \\frac{1}{L_t} M_0(dx_0)\\prod_{s=1}^t M_s(x_{s-1},dx_s)\n", "\\prod_{s=0}^t G_s(x_{s-1}, x_s)\n", "\\end{equation}\n", "\n", "where: \n", "\n", "* $M_0(dx_0)$ is the $N(0,1)$ distribution; \n", "* $M_s(x_{s-1},dx_s)$ is the $N(x_{s-1}, 1)$ distribution; \n", "* $G_s(x_{s-1}, x_s)= \\mathbb{1}_{[0,\\epsilon]}(x_s)$\n", "\n", "Let's define the corresponding `FeymanKac` object:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "import seaborn as sb\n", "import numpy as np\n", "from scipy import stats\n", "\n", "import particles\n", "\n", "class GaussianProb(particles.FeynmanKac):\n", " def __init__(self, a=0., b=1., T=10):\n", " self.a, self.b, self.T = a, b, T\n", "\n", " def M0(self, N):\n", " return stats.norm.rvs(size=N)\n", "\n", " def M(self, t, xp):\n", " return stats.norm.rvs(loc=xp, size=xp.shape)\n", " \n", " def logG(self, t, xp, x):\n", " return np.where((x < self.b) & (x > self.a), 0., -np.inf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The class above defines: \n", "\n", "* the initial distribution, $M_0(dx_0)$ and the kernels $M_t(x_{t-1}, dx_t)$, through methods `M0(self, N)` and `M(self, t, xp)`. In fact, these methods simulate $N$ random variables from the corresponding distributions. \n", "* Function `logG(self, t, xp, x)` returns the log of function $G_t(x_{t-1}, x_t)$.\n", "\n", "Methods `M0` and `M` also define implicitly how the $N$ particles should be represented internally: as a (N,) numpy array. Indeed, at time $0$, method `M0` generates a (N,) numpy array, and at times $t\\geq 1$, method `M` takes as an input (`xp`) and returns as an output arrays of shape (N,). We could use another type of object to represent our $N$ particles; for instance, the `smc_samplers` module defines a `ThetaParticles` class for storing $N$ particles representing $N$ parameter values (and associated information).\n", "\n", "Now let's run the corresponding SMC algorithm:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEMCAYAAAA1VZrrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXRURfr/8XdVOhJ2yAIRhOEHBDGiICYiGQSViA4yiMo+qLgATthkkUUkASIYlrBv+pUdZBHI4MKMGhdQQQGVcUFGgRFFAiEkEFGQJFW/PzJmZIhJZ73dned1jud0N327P8+5wtO37r1VylprEUIIIQqgnQ4ghBDC80mzEEIIUShpFkIIIQolzUIIIUShpFkIIYQolDQLIYQQhXI5HcAd+/fvZ8WKFRhj6NixI926dXM6khBCVCgef2RhjGHZsmU89dRTzJkzhw8++IBjx445HUsIISoUj28Whw4dIjQ0lLp16+JyuYiKimLv3r1OxxJCiArF44eh0tPTCQoKynseFBTEN998U+A2x48fL/b3BQcHk5aWVuztPY3U4/l8rSZfqwd8r6b86qlXr16B23h8s8hvNhKl1CXPk5OTSU5OBiAhIYHg4OBif5/L5SrR9p5G6vF8vlaTr9UDvldTcerx+GYRFBTE6dOn856fPn2a2rVrX/Ke6OhooqOj856X5BdARfgF4c18rR7wvZp8rR7wvZqKc2Th8ecsmjRpQkpKCqmpqWRnZ7Nr1y4iIiKcjiWEEBWKxx9Z+Pn58cgjjzB16lSMMdx22200aNDA6VhCCFGheHyzAGjdujWtW7d2OoYQQlRYHj8MJYQQwnnSLIQQQhRKmsV/2KyLmA3/R86pE05HEUIIjyPN4lf//ga783XShvTGJK3FXvjZ6URCCOExpFn8h2p2LTp+CQE334rdvgnz9F8x772BNTlORxNCCMdJs/gNFRRCzRGT0ONnQnBd7OqFmPgR2K/+6XQ0IYRwlDSLfKjGV6PHTkcPGgPnf8bMnkjOwmewJ2S2WyFExSTN4ncopVAR7dDxi1H3PwT/+hwzaShm/fPYc5lOxxNCiHLlFTflOUn5X4G6635sVEfsyy9i39mO/fAdVNe+qNs6o7Sf0xGFEKLMyZGFm1SNWuh+Mei4+dCoGXbD/2GmPYk9etjpaEIIUeakWRSRqt8Q/cQk1MAxkJGGmToKs2kZ9sJ5p6MJIUSZkWGoYlBKoSLbYa9thd2yGvvmNuzHu9B9H0e1jHQ6nhBClDo5sigBVaUa+oEY9NgEqBSAWRhPztIE7Jl0p6MJIUSpkmZRClTTcHTsXFS3fvDPvZjYGMy727HGOB1NCCFKhTSLUqJc/ui7e6InLYA/NMWuW4qZMQ77w1GnowkhRIlJsyhlqm499Mh41MNPwMkfMPFPYLauxl78xeloQghRbHKCuwwopVBRt2Ovi8BuXoH9+2bsvvfR/f6KCr/B6XhCCFFkcmRRhlT1GuiHh6NHPQNKY+bEYV5IxGaecTqaEEIUiTSLcqCaX4+eNB/VpRd23weY2MGY99/EWut0NCGEcIs0i3Ki/K9A3/MXdOxcuLIBdtUCzKwJMjmhEMIrSLMoZ6peQ/ST01APDoFj/8ZMHoZ5eT02K8vpaEII8bukWThAaY2+pVPujLato7CvrMdMGYb9+kunowkhRL6kWThI1aiNHjAaPTwOsrMxs57C/G0tNjvb6WhCCHEJaRYeQLW4ER03HxXVEfvaJszM8dhTJ5yOJYQQeTz6PotNmzbx1ltvUaNGDQD69OlD69atHU5VNlRAZVT/YZhrW2PXLMJMGY7qF4Nu08HpaEII4dnNAuDuu++ma9euTscoNzqyHbZxs9z7MV5IxHzxCarvIFTlKk5HE0JUYDIM5YFUUB306GmoP/fBfrQDE/8E9t9fOx1LCFGBeXyzeP311xk9ejSLFy/m3LlzTscpN8rPD921D3rMNMjJwUwfi9n+EtbkOB1NCFEBKevwbcTx8fGcOXP59Be9e/cmLCws73zFxo0bycjIICYm5rL3Jicnk5ycDEBCQgIXL14sdh6Xy0W2h12NZH76kcwlM/jlg7fwb9GamsNj8Quu49a2nlhPSfhaPeB7NflaPeB7NeVXzxVXXFHgNo43C3elpqYyffp0EhMTC33v8ePHi/09wcHBpKWlFXv7smKtxe56G7v+OXD5o3o8goq6HaVUgdt5aj3F5Wv1gO/V5Gv1gO/VlF899erVK3Abjx6GysjIyHu8Z88eGjRo4GAaZyml0H/siJ44F0LrY1fOw8wYL+tlCCHKhUdfDbV27Vq+/fZblFKEhIQwcOBApyM5TtWthx6TgP0gGbtlFSb+CVR0V1SX3qiAyk7HE0L4KI9uFkOHDnU6gkdSWqNu6YRtdTN26yrs60nYve+hew2AG24udGhKCCGKyqOHoUTBVPUa6IeGoscmQOWqmCXPYhbEy93fQohSJ83CB6im4ein56B6PAJff4mJG4J5bZPMZCuEKDXSLHyEcrnQnbqhpyyC6yOwf1uLmTKMXz7b53Q0IYQPkGbhY1RgMH6Pj0MPi4OcHM7EDcMsnY5NP+V0NCGEF/PoE9yi+NR1N6KbL6Tye6/z0+ZV2M/3oe7uibqjG8rf3+l4QggvI0cWPkz5X0G1ng/nDk1dewM2aQ1m0lDsFx87HU0I4WWkWVQAKrgufjFPoYdPAsDMm0zOoqly1ZQQwm3SLCoQ1aI1etIC1H0PwVf/zL1q6uUXsRd/cTqaEMLDSbOoYJS/P/pP96OnLEa1aoN9ZQMmdjB2/4d4yTRhQggHSLOooFRgMHrgk+iR8XBFJcyiaZhZE7CHDzodTQjhgaRZVHDqmpbo2HmoPgMh5XtMwhhyFk3DpnzvdDQhhAeRS2cFyuVC3d4FG9URm7wN+3oSJm4PKuo2VNe+qMAQpyMKIRwmzULkUQGVUV16Yzt0xm5/Cfvua9iPdqJuvxv1p+6oajWcjiiEcIgMQ4nLqOo10L0eRT/zHOqm9tg3X8Y8NTB3vqlfLjgdTwjhAGkW4nepoBD0w8PRcfPh6uty55t6aiDmvTfkyikhKhhpFqJQqn5D/AZPQI+bAXXqYVcvxMyJxaaddDqaEKKcSLMQblNNmqPHPIt6IAaOfI2ZNBTz7nasMU5HE0KUMWkWokiUUuj2d6EnL4AmzbHrlmJmT5SpQ4TwcdIsRLGooDroJyajHhwCRw9hJg/DvPOaHGUI4aOkWYhiU0qhb+mEnrQQml6DffE5TOLT2NQUp6MJIUqZNAtRYiooBD18EuqhofD9kdyjjLdekaMMIXyINAtRKpRS6HZ35B5lNGuB3fB/mFlPYVOOOR1NCFEKpFmIUqUCg9HDYlH9h8MPR3OPMpLWYH+RadCF8GbSLESpU0qh/9gRHb8YddMt2O0vYeIGY/+51+loQohikmYhyoyqURv9yAj06Gm506AvjM+d0fb0KaejCSGKyPFmsXv3bkaOHEmvXr04fPjwJX+WlJTE0KFDGT58OPv373cooSgpdXULdOzc3BX6DnyKiY3B/GMLNjvb6WhCCDc53iwaNGjA6NGjueaaay55/dixY+zatYvZs2czYcIEli1bhpGra7yWcv26Qt8iCG+F3bIKM2U49usvnI4mhHCD483iqquuol69epe9vnfvXqKiovD396dOnTqEhoZy6NAhBxKK0qSC6uTOMzXkabj4C2bmU5jlc7GZZ5yOJoQogMeuZ5Genk5YWFje88DAQNLT0/N9b3JyMsnJyQAkJCQQHBxc7O91uVwl2t7TeGw9HTtj293OuZdW8vO2F+HzvVTt91cq39EVpX//N4zH1lMCvlaTr9UDvldTcepxu1ls376ddu3aUaNG0RfAiY+P58yZy3859u7dm8jIyHy3KcoU2NHR0URHR+c9T0tLK3LGXwUHB5doe0/j8fXc1R3dsg1m3VJ+XDqDH1//G7rfX1ENm+T7do+vpxh8rSZfqwd8r6b86slvhOe33G4Wn3/+OevXr+faa6+lffv2REZG4u/v79a2EydOdPdr8gQFBXH69Om85+np6QQGBhb5c4TnU1c2QI96BvvRDuymZZhnRuWuznfPX1CVqzgdTwhBEc5ZjB07lsWLF9OqVStee+01Bg4cyNKlSzlw4ECZBIuIiGDXrl1kZWWRmppKSkoKTZs2LZPvEs5TSqFvvhX9zBLUrXdh334VMzEGs2enLLQkhAdQtph/E48ePcrChQv57rvvCA4OpmPHjnTu3JmAgIAifc6ePXtYvnw5mZmZVK1alUaNGjFhwgQAtm7dyjvvvIPWmv79+3PDDTe49ZnHjx8vcj2/qgiHm97AfvsNZu0SOHoIrmmJ7vs4KrS+19ZTEF+rydfqAd+rqTjDUEVuFp9//jnvvfcee/fupUmTJnTo0IHg4GC2b9/O2bNnmTJlStGTlzJpFv/lzfVYk4Pd8To2aQ1k/YK6635C+j3O6R9/dDpaqfLmfZQfX6sHfK+mMj1nsXr1anbt2kWVKlVo3749iYmJl5xDCAsL4+GHHy5iZCF+n9J+qNs6Y1u3xW5egX11I6f3voe990Fo3RallNMRhagw3G4WWVlZjB49+nfPG7hcLhISEkotmBC/UjVrox4dif1jNGxahlmaAE2ao7v3RzUNdzqeEBWC2ye4lVL5NoqVK1fmPa5fv36phBIiP6r59QTNWZ27Ol9aKmb6OHIWT8OekGnQhShrbjeLHTt25Pv6zp07Sy2MEIVRfn65q/NNXYq65y9w4J+YuCGYdUuxmRlOxxPCZxU6DPX2228DkJOTk/f4V6mpqVSvXr1skglRAFUpANWlF7b9ndhXN2B3vo7d/Q7qrntRd3RDVSraVXlCiIIV2izee+89ALKzs/Me/6pmzZoMHjy4bJIJ4QZVoxaq7+PY2/+MSVqN3fYi9t1/oO7pi4rqiPLzczqiED6h0GYRFxcHwIYNG+jdu3eZBxKiOFRoffz+Oh576ABm80rs6oXYd7ejHxmBqv8Hp+MJ4fUKPGfx21swevbsiTEm3/+E8BSqaTh67HT0oDGQcRrzzIjctTNMjtPRhPBqBR5Z9O/fn1WrVgHQp0+f333fxo0bSzeVECWglIKIduirr8OsXYzdsgq7/yP0I0+g6hR845EQIn8FNovExMS8xwsXLizzMEKUJlW9JvrxcbkTFK5/DjN5OKp7f1SHPxU4DboQ4nIFNovfznceEhJS5mGEKG1KKdTNt2Kvvg6zaj72xeewn36I7j8MFSj/TwvhrgKbxYIFC9yaUmHIkCGlFkiIsqBqB6GHT8K+9zp203LMpKGo3gNQbW+XaUOEcEOBzSI0NLS8cghR5pRSqPZ3Ya9phVkxF7tiXu5RxgMxqBq1nY4nhEcrsFn06NGjvHIIUW5USCh69FRs8svYpLWYuKG5Cy21uwPl8tiVhoVwVIF/Mw4cOEB4eO5EbV988cXvvq9Fixalm0qIMqa0H6rTvdgWN+ZeMbVuCfatl9H3PQSt2sjQlBD/o8BmsWzZsrwropYsWZLve5RScqWU8FqqXkP0k8/CP/dgtqzCLJ4GTcNzZ7Rt0tzpeEJ4jGKvlOfJZPGj/5J63GdzcrDvv4l9+UXIPAM3RqHve7DM782QfeT5fK2mMl38CMAYw9dff01GRgaBgYGEhYWh5Xp14SOUnx+qw13YNh2wb/wN+0YSZv9HufdldOmFql7T6YhCOMbtZnH06FFmzpxJVlYWgYGBpKen4+/vz+jRo2nUqFEZRhSifKmAyqiufXJntH1lA/bd7dhdb6Huuh8VfQ+qUiWnIwpR7txuFkuWLOHOO++kS5cuKKWw1vLaa6+xZMkSpk+fXpYZhXCEqhWIeiAGG/1nzJZV2L+txb73Ru4Nfc2vdzqeEOXK7TGklJQU7r777ryrRJRSdO7cmRMnTpRZOCE8gbqyAX5DnkaPngp+fpjEpzHrlmAvnHc6mhDlxu1mccMNN7Bv375LXtu3bx833HBDqYcSwhOpq69Dx85HRd+D3fEPzKSh2K/+6XQsIcqF29N9GGOYO3cujRs3JigoiNOnT3PkyBEiIiLKJagQnkBVqoTq9Sj2xraYFfMxsyei2t+F6tEfFVDF6XhClJkiTffRoEGDvMdXXXUVLVu2LJtUQng41TQcHTcPu20d9s1t2C8+Rj80FBXeyuloQpQJx6f72L17Ny+99BI//PAD06ZNo0mTJkDu+t4jRozIu/Y3LCyMgQMHlnkeIdylrqiE6vEItnUUZuU8zJxY1C2dUD0eQVWWowzhW4p0n0V2djbHjx8nMzPzktdLMt1HgwYNGD16NM8///xlfxYaGsrMmTOL/dlClAfVpDl64lzsyy9i39iG/fIT9ANDUC1aOx1NiFLjdrM4ePAgs2fPJisri/Pnz1O5cmUuXLhAUFBQiab7uOqqq4q9rRCeQl1RCdX94f8cZczHzJuUOzFhj0dQVao6HU+IEnP7aqhVq1bRtWtXVqxYQeXKlVmxYgX3338/nTp1KrNwqampjBkzhri4OL766qsy+x4hSotqfDV64hzUXfdjP3gr94qpzz92OpYQJeb2kcXx48fp3LnzJa9169aNwYMH07Vr1wK3jY+P58yZM5e93rt3byIjI/Pdpnbt2ixevJjq1atz5MgRZs6cSWJiIlWqXD4WnJycTHJyMgAJCQmXrPBXVC6Xq0TbexqpxyGDRpF12584u3AqOfMnE3D73VR/ZBi6avXL3uo1NbnJ1+oB36upOPW43SyqVKnC+fPnqVq1KrVq1eLYsWNUq1aNCxcuFLrtxIkTixQKwN/fH39/fwAaN25M3bp1SUlJyTsB/lvR0dFER0fnPS/JhF8VYcIwb+ZV9QTWwY6fhXp1Axf+sYULn+xGPzAYdf2lP5C8qiY3+Fo94Hs1FWciQbeHodq0acOnn34KwO23387kyZMZN24cbdu2LUbUwmVmZmKMAeDkyZOkpKRQt27dMvkuIcqK8vdH3/sA+qlZULU6ZkE8Zvlc7E/nnI4mRJEUe4rygwcPcv78eVq2bFmimWf37NnD8uXLyczMpGrVqjRq1IgJEybw4YcfsmnTJvz8/NBa06NHD7dvAJQpyv9L6vEcNisLu30TdvtLUL1W7nKuLW/y6pry42v1gO/VVJwjiyI3i/T0dNLT0wkMDCQwMLDoKcuBNIv/kno8jz16GLNyHhz7FnXzrQTHjCP9l4tOxyo1vrCP/pev1VSm61mkpaUxf/58vv76a6pVq8a5c+cICwtj2LBhhISEFC+xEBWQ+kMT9IRE7PaXsNtfIn3MN9gnJqOC6jgdTYjf5fb40aJFi2jcuDErV67khRdeYOXKlTRp0oRFixaVZT4hfJJy+aO79kU/+Szmx0zMjHHYk8U/IhairLndLI4cOUK/fv0ICAgAICAggH79+nHkyJEyCyeEr1NNmlM7fgFcvIiZOR77w3dORxIiX243i7CwMA4dOnTJa4cPH6ZZs2alHkqIisT//zVDPzkNUJhZ47FHDzsdSYjLFHjOYuPGjXmP69aty7PPPkvr1q3zpij/9NNPadeuXZmHFMLXqXoN0WOmYRInYhKfRg+PQzVp7nQsIfIUeGRx+vTpvP+ysrJo06YN/v7+ZGZm4u/vz0033cTFi75zFYcQTlJ16qHHJED1Gpg5sdh/fe50JCHyFHhkERMTU145hBCACgrJPek9eyJm3mR0zHhUixudjiWE++csIHcd7s2bN/P888+zefNmUlJSyiqXEBWWqhWIfvJZuPIqzMKp2E92Ox1JCPebxb59+xg3bhw//PAD1apV4/jx44wbN+6ydbmFECWnqtdAj3oGGjXFPDcd89EOpyOJCs7tm/LWr1/Pk08+eclCR19++SXLly+XdbiFKAOqSjX0E5MxC5/BLpuNufgL+payWxJAiIK4fWSRnp7ONddcc8lrzZs35/Tp06UeSgiRSwVURg+LhWtbY1cvxKxeiP1ZJiEU5c/tZtGoUSNeeeWVS1579dVXadSoUWlnEkL8hrqiEnrwU6g778W+n4yJHYL99EOnY4kKxu1hqMcee4zp06fz97//Pe8+i0qVKjFmzJiyzCeEIHd6ENX9YWzkLZiVCzCLp8GNUeg+g1A1azsdT1QAbjULYwwZGRnMmDGDb7/9Nm/W2aZNm+Jyud1vhBAlpP7QNHcSwjeSsK9swHz1Garno6io21FKOR1P+DC3hqG01syYMYOAgACaN29OVFQUzZs3l0YhhAOUy4Xu3AMdNw/qNcSunIeZG4c9dcLpaMKHuX3O4pprruHrr78uyyxCiCJQoVehn5yG6vs4HP4XZtJQzJvbsCbH6WjCB7l9aBASEsKzzz5LREQEQUFBlxzy9urVq0zCCSEKprRG3dYZe30kZt0S7KZl2L3voR8aiqr/B6fjCR/i9pHFxYsXiYyMRClFenr6JfNGCSGcpYJC0EMnoh4bBadOYOJHYLa9iM3Kcjqa8BFuH1nIPFFCeDalFKpNB2x4K+zGF7CvbsB+/EHuUYbMYCtKqEhnqFNSUti9e3fe1VBt27blyiuvLKtsQohiUNVroh4bhW3TAbN2MWb6WNTtXVDd+qECKjsdT3gpt4eh3n//fcaMGcPRo0cJCAjgu+++Y+zYsbz//vtlmU8IUUzqugj05IWoW/+EfesVzKSh2C8/dTqW8FJuH1ls2LCB8ePHEx4envfaV199xcKFC2UBJCE8lAqogur7ODayPWb1AszcOFTb21G9HkVVre50POFF3D6yOH/+/GVLqIaFhXHhwoVSDyWEKF0qLBwdOw/VuSd2zw7MxBjsvvex1jodTXgJt5tFly5dWL9+fd7KeBcvXmTDhg106dKlzMIJIUqP8r8CfW8/9ITZEBiCeW4GZvE07OlTTkcTXsDtYag33niDM2fOsH37dqpVq8a5c7kzX9aqVYs33ngj731LliwpUoA1a9bw8ccf43K5qFu3LjExMVStWhWApKQk3n77bbTWPPzww7Rq1apIny2EuJxq8P/Q42dik1/GvrwOE/tXVOeeqE73ovz9nY4nPJTbzWLo0KFlEuD666+nb9+++Pn5sXbtWpKSkujXrx/Hjh1j165dzJ49m4yMDOLj45k3bx5aF2lxPyFEPpSfX+4sthF/xGxajv3bWuyut9C9B6Cuk/VpxOXcbha/PbFdmlq2bJn3uFmzZnz4Ye7Uy3v37iUqKgp/f3/q1KlDaGgohw4duuy8iRCi+FRQHfz+Og574FPM+ucx86dAy5vQvR5DhYQ6HU94kGL9TH/ooYdKOwcAb7/9dt5QU3p6OkFBQXl/FhgYSHp6epl8rxAVnQq/AR03H9W9Pxz8DBM7OPcO8F9+cTqa8BDFmja2qFdQxMfHc+bMmcte7927N5GRkQBs3boVPz8/brnlliJ/R3JyMsnJyQAkJCQQHBxcpHy/5XK5SrS9p5F6PJ9H1fSXgeTcdS/nVi3kwqsbUB+9S7VHh1PppvZuT4HuUfWUEl+rqTj1lMsc4xMnTizwz999910+/vhjYmNj8/6H/HWBpV/9etd4fqKjo4mOjs57npaWVuyswcHBJdre00g9ns/zalLw4FB0m9sw65/jbMJ4uPaG3PMZoVcVurXn1VNyvlZTfvXUq1evwG2KNQw1e/bs4myWr/3797Nt2zbGjh1LpUqV8l6PiIhg165dZGVlkZqaSkpKCk2bNi217xVCFExd3QL99BxUr8fgyL8wcUMw65ZiMzOcjiYc4PaRxcmTJ/N97u/vT61atYp9ldKyZcvIzs4mPj4eyL3Rb+DAgTRo0IC2bdsycuRItNY8+uijciWUEOVMuVyo6K7Ym9pjX9mA3fkP7O53UHfdi7qjG6pSgNMRRTlR1s2TAwWtWaG15sYbb+Sxxx6jVq1apRauuI4fP17sbSvC4aY387V6wLtqsieOYZLWwCe7oWYg6p6+qKiOKD+/vPd4Uz3u8rWaijMM5faRxaBBgzhw4ADdu3fP+6LNmzdz9dVXEx4ezrp161i2bBmjRo0qXnohhMdToVfh99fx2EMHMJtXYlcvxCa/jL7/IbguQtYB92Fuj+ts2rSJgQMHEhoaisvlIjQ0lAEDBrBlyxbq169PTEwMBw4cKMusQggPoZqGo8dORz8+DrKzMAviMYlPY7/9xuloooy43SystZw6dekcMmlpaRhjAAgICCAnR9b+FaKiUEqhboxCT16E6jMQfjiKmTqKs3MmYX/MdDqeKGVuD0N17tyZKVOmcOuttxIUFER6ejrvvPMOnTt3BuCTTz6Ru6uFqICUy4W6vQu27e3Yf2zlwhtJ8Pkn6EFjZIU+H+L2CW7Ivcx19+7dZGRkUKtWLaKiojxycj85wf1fUo/n87Waap49TXrCOMhIQ3Xvj+rY1evPZfjaPirTE9wArVq18sjmIITwHP5NrkZPnINZMR+7cRn2mwPoh4ahqlR1OpooAbebRXZ2Nlu3bmXnzp1kZGRQu3Zt2rdvz3333YfLVS43ggshvISqUg0dMx775jbs1lWYZ0agHx+LatjE6WiimNz+V37t2rUcPnyYAQMGEBISwqlTp9iyZQs///wz/fv3L8OIQghvpJRCdeqGbdwM89xMzLNjUH0GoG650+uHpSoit6+G+vDDDxkzZgwtW7akXr16tGzZktGjR7N79+6yzCeE8HKqaTg6di40a4Fdsxi7fC72F1mO2dsU6dJZIYQoDlW9Jnp4LKprX+xH72KmjsKmfO90LFEEbjeLtm3bMn36dPbv38+xY8fYv38/M2fO5Oabby7LfEIIH6G0H/rPvdFPTIZzmZipozAfJMsPUS/h9jmLfv36sWXLFpYtW0ZGRgaBgYFERUXRvXv3sswnhPAxKrwVeuJczAuzsCvnYz98F90vBlW34Es3hbMKbBZffPHFJc+vvfZarr32Wqy1eSeoDh48SIsWLcouoRDC56jaQehRU7E7X8duXY2ZNBR1d0/UXfehXP5OxxP5KLBZLFmyJN/Xf20UvzaNhQsXln4yIYRPU1qjbv0TtlUb7MYXsNvWYffsRD8wGBUW7nQ88T8KbBaLFi0qrxxCiApK1QpEDRqDjbods24pZsY4VPs7Ufc9hKpazel44j9kNSEhhEdQ10WgJy/MvTfjvTcxsTGYPTvlBKD9zlgAAA7QSURBVLiHkGYhhPAYqlIAuscj6KcToXYw9v9mYeZPxp464XS0Ck+ahRDC46iGTdBPzcxd//ubA5hJQzHvbpejDAdJsxBCeCSl/dDRXdFTFkHTcOy6pZgF8dizGU5Hq5CkWQghPJoKDEEPj0P1HggHP8NMGord/6HTsSocaRZCCI+ntEZ37IJ+ejbUDsIsmoZZtQB74bzT0SoMaRZCCK+h6jVEPzUL9af7sR8kY6YMxx4+6HSsCkGahRDCqyiXP/q+h9Cjp4IxmOnjMNvWYbOznY7m06RZCCG8kmrWAh07D3VzB+yrGzHTx2JP/OB0LJ/l+BJ3a9as4eOPP8blclG3bl1iYmKoWrUqqampjBgxIm9d2LCwMAYOHOhwWiGEJ1FVqqIeGYG9PhKzZjEm/glU7wHoWzo5Hc3nON4srr/+evr27Yufnx9r164lKSmJfv36ARAaGsrMmTMdTiiE8HQqoh26yTWYFXOxqxdijvwL1XcQyv8Kp6P5DMeHoVq2bImfnx8AzZo1Iz093eFEQghvpGoHoZ+YhOrcE/v+m5jp47BpJ52O5TMcbxa/9fbbb9OqVau856mpqYwZM4a4uDi++uorB5MJIbyB0n7oe/uhB0+A1OOYZ0Ziv/zU6Vg+QdlyuH8+Pj6eM2fOXPZ67969iYyMBGDr1q0cPnyY0aNHo5QiKyuLCxcuUL16dY4cOcLMmTNJTEykSpUql31OcnIyycnJACQkJHDx4sViZ3W5XGT70FUVUo/n87WaPKWe7OPfc3b6eLK//zdV+wyg6v0PonTxfh97Sk2lJb96rrii4CG7cmkWhXn33Xd58803iY2NpVKlSvm+Z9KkSTzwwAM0adKk0M87fvx4sbMEBweTlpZW7O09jdTj+XytJk+qx/5yAbt6EXbPDmh5E/qRJ1BVij7tuSfVVBryq+fXi4l+j+PDUPv372fbtm2MHTv2kkaRmZmJMQaAkydPkpKSQt26dZ2KKYTwQqpSAOqxkblThXzxMWbqKOyxb52O5ZUcvxpq2bJlZGdnEx8fD/z3EtkDBw6wadMm/Pz80FozYMAAqlWThVCEEEWjlEJ17IJt2Bjz3HTMs0+iHhyCbtPB6WhexSOGoUqbDEP9l9Tj+XytJk+ux55Jxzw/A745gOr4Z9T9/VH+ha/57ck1FYdXDkMJIUR5UbUC0SOfQUV3xb71CubZ0dgfvnM6lleQZiGEqFCUy4Xu9Vju5bVn0jHPjMAkv4z9zzlSkT9pFkKICkm1aoOeNB/CW2E3voCZNwmbcdrpWB5LmoUQosJSNWqjhzyN6hcDh77KXVhp3/tOx/JI0iyEEBWaUgrd4S70xLlQtx7muRmY5XOwP//kdDSPIs1CCCEAFVofPSYB1aU39sMduQsrff2l07E8hjQLIYT4D+Vyoe/pix6bAFpjZj2F2boKm5XldDTHSbMQQoj/oZo0z11Yqd0d2L9vIX3sYxX+zm9pFkIIkQ8VUBn94BD04AmY9DTM1JGYv2/BmhynozlCmoUQQhRAtWpD0Px10PIm7NZVmBnjsSeLP0uEt5JmIYQQhdA1aqEHjUU9NgpSvsdMGY5557UKdSOfNAshhHCDUgrdpgN60kJodi32xecwc+Owp085Ha1cSLMQQogiULWD0MPiUA/EwJF/YSYPxXzwFj44J+slpFkIIUQRKaXQ7e9Cx82HBv8Pu3IeZtFUbGaG09HKjDQLIYQoJhUSih41FdXzUfjyU0zcEMzud3zyKEOahRBClIDSGn3HPejYuVC3Pnb5HMycWGyqb10xJc1CCCFKgbqyQe50IX95HL79BhM3FPPaJmy2b9z9Lc1CCCFKidIafWtn9JRFqJY3Yf+2FhM/AnvogNPRSkyahRBClDJVKwj9+Fj0kIlw4Txm+jjMmsXYn885Ha3YXE4HEEIIX6VaRqKvboHd9iL2rVew//wI1WsAKuKPKKWcjlckcmQhhBBlSAVURvd6FD1hFtQMxD4/A7MgHnvqhNPRikSahRBClAP1h6bop2ahej0KX3+BiR2cO/35+Z+djuYWaRZCCFFOlJ8fOvoe9JTFqIh22L9vwUwYhNnxD2yOZ89mK81CCCHKmQoMRj86Aj0hEULrY9cuxsQ/gf3yU6ej/S7HT3Bv2LCBffv2oZSiZs2axMTEEBgYCEBSUhJvv/02WmsefvhhWrVq5XBaIYQoPapRGPrJZ+HT3ZjNKzFz46DFjegeD6PqNXQ63iUcP7Lo2rUrs2bNYubMmbRu3ZrNmzcDcOzYMXbt2sXs2bOZMGECy5Ytw1Sg6YCFEBWDUgrVOgo9eRGqx8Nw+CBm8jDMuiXYH886HS+P482iSpUqeY9/+eWXvMvJ9u7dS1RUFP7+/tSpU4fQ0FAOHTrkVEwhhChTyt8f3ele9NTnUB3+hN35eu75jNeTPOJ8huPDUADr169n586dVKlShbi4OADS09MJCwvLe09gYCDp6elORRRCiHKhqtdA9R2Eva0zZvNK7OYV2E93ox8bhQqu61wuWw7TI8bHx3PmzJnLXu/duzeRkZF5z5OSksjKyqJnz5688MILNGvWjPbt2wOwZMkSbrjhBm6++ebLPic5OZnk5GQAEhISuHjxYrGzulwusrOzi729p5F6PJ+v1eRr9YCzNV14700yl84AoMbjYwi45Y4Sf2Z+9VxxxRUFblMuzcJdp06dIiEhgcTERJKSkgC49957AZg6dSo9evSgWbNmhX7O8ePFn+0xODiYtLS0Ym/vaaQez+drNflaPeB8TTbtJOaFRDh8EBXVEdVnACqgSuEb/o786qlXr16B2zh+ziIlJSXv8b59+/ICR0REsGvXLrKyskhNTSUlJYWmTZs6FVMIIRyjguuin3wW1aU3dvc7uZMT/vubcs3g+DmLdevWkZKSglKK4OBgBg4cCECDBg1o27YtI0eORGvNo48+itaO9zYhhHCE8vND3dMXe01LzLJEzPQxqG79UJ3uRZXDv40eNQxVWmQY6r+kHs/nazX5Wj3geTXZn85h1iyEj3fBNS3RjzyBqhXk9vZeOQwlhBCiaFTVauhBY1EPDsm7L8Pu/6hMv1OahRBCeCGlFPqWTuin50BgCGbRVMxLy8vs+6RZCCGEF1NXXoUeNxPVqRuEhJbZ9zh+glsIIUTJKH9/VI9HyvQ75MhCCCFEoaRZCCGEKJQ0CyGEEIWSZiGEEKJQ0iyEEEIUSpqFEEKIQkmzEEIIUShpFkIIIQrlkxMJCiGEKF1yZPE/xo0b53SEUiX1eD5fq8nX6gHfq6k49UizEEIIUShpFkIIIQrlN2nSpElOh/A0jRs3djpCqZJ6PJ+v1eRr9YDv1VTUeuQEtxBCiELJMJQQQohCyXoW/7F//35WrFiBMYaOHTvSrVs3pyOV2ODBgwkICEBrjZ+fHwkJCU5HKpLFixfzySefULNmTRITEwE4d+4cc+bM4dSpU4SEhDBixAiqVavmcFL35VfTpk2beOutt6hRowYAffr0oXXr1k7GdFtaWhqLFi3izJkzKKWIjo6mc+fOXruffq8eb95HFy9eJC4ujuzsbHJycrj55pvp2bNn0feRFTYnJ8cOGTLEnjhxwmZlZdnRo0fb77//3ulYJRYTE2PPnj3rdIxi+/LLL+3hw4ftyJEj815bs2aNTUpKstZam5SUZNesWeNUvGLJr6aNGzfabdu2OZiq+NLT0+3hw4ettdb+/PPPdtiwYfb777/32v30e/V48z4yxtjz589ba63Nysqy48ePt//617+KvI9kGAo4dOgQoaGh1K1bF5fLRVRUFHv37nU6VoUXHh5+2S+dvXv30qFDBwA6dOjgdfspv5q8We3atfNOlFauXJn69euTnp7utfvp9+rxZkopAgICAMjJySEnJwelVJH3kQxDAenp6QQFBeU9DwoK4ptvvnEwUemZOnUqAHfccQfR0dEOpym5s2fPUrt2bSD3L3ZmZqbDiUrH66+/zs6dO2ncuDEPPvigVzaU1NRU/v3vf9O0aVOf2E+/refgwYNevY+MMYwdO5YTJ05w5513EhYWVuR9JM0CsPlcEKaUciBJ6YqPjycwMJCzZ8/yzDPPUK9ePcLDw52OJf5Hp06d6N69OwAbN25k9erVxMTEOJyqaC5cuEBiYiL9+/enSpUqTscpsf+tx9v3kdaamTNn8tNPPzFr1iy+++67on9GGeTyOkFBQZw+fTrv+enTp/M6rjcLDAwEoGbNmkRGRnLo0CGHE5VczZo1ycjIACAjIyPvhKM3q1WrFlprtNZ07NiRw4cPOx2pSLKzs0lMTOSWW26hTZs2gHfvp/zq8fZ99KuqVasSHh7O/v37i7yPpFkATZo0ISUlhdTUVLKzs9m1axcRERFOxyqRCxcucP78+bzHn332GQ0bNnQ4VclFRESwY8cOAHbs2EFkZKTDiUru17+wAHv27KFBgwYOpikaay1Lly6lfv36dOnSJe91b91Pv1ePN++jzMxMfvrpJyD3yqjPP/+c+vXrF3kfyU15//HJJ5+watUqjDHcdttt3HfffU5HKpGTJ08ya9YsIPekVrt27byuprlz53LgwAF+/PFHatasSc+ePYmMjGTOnDmkpaURHBzMyJEjvWrsOL+avvzyS7799luUUoSEhDBw4ECvObI9ePAgsbGxNGzYMG/otk+fPoSFhXnlfvq9ej744AOv3UdHjx5l0aJFGGOw1tK2bVu6d+/Ojz/+WKR9JM1CCCFEoWQYSgghRKGkWQghhCiUNAshhBCFkmYhhBCiUNIshBBCFEqahRBCiEJJsxCinAwePJjPPvvM6RhCFIs0CyGEEIWSm/KEKAcLFizg/fffx+VyobWme/fu3HPPPU7HEsJt0iyEKCeDBw9m0KBBXH/99U5HEaLIZBhKCCFEoaRZCCGEKJQ0CyGEEIWSZiFEOalVqxapqalOxxCiWOQEtxDlZO/evSxfvpzz589z33330bVrV6cjCeE2aRZCCCEKJcNQQgghCiXNQgghRKGkWQghhCiUNAshhBCFkmYhhBCiUNIshBBCFEqahRBCiEJJsxBCCFEoaRZCCCEK9f8Bewqi/8FtqXAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fk_gp = GaussianProb(a=0., b=1., T=30)\n", "alg = particles.SMC(fk=fk_gp, N=100)\n", "alg.run()\n", "\n", "plt.style.use('ggplot')\n", "plt.plot(alg.summaries.logLts)\n", "plt.xlabel('t')\n", "plt.ylabel(r'log-probability');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That was not so hard. However our implementation suffers from several limitations: \n", "\n", "1. The SMC sampler we ran may be quite inefficient when interval $[a,b]$ is small; in that case many particles should get a zero weight at each iteration. \n", "2. We cannot currently run the SQMC algorithm (the quasi Monte Carlo version of SMC); to do so, we need to specify the Markov kernels $M_t$ in a different way: not as simulators, but as deterministic functions that take as inputs uniform variates (see below). \n", " \n", "Let's address the second point: " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "class GaussianProb(particles.FeynmanKac):\n", " du = 1 # dimension of uniform variates\n", " \n", " def __init__(self, a=0., b=1., T=10):\n", " self.a, self.b, self.T = a, b, T\n", "\n", " def M0(self, N):\n", " return stats.norm.rvs(size=N)\n", "\n", " def M(self, t, xp):\n", " return stats.norm.rvs(loc=xp, size=xp.shape)\n", " \n", " def Gamma0(self, u):\n", " return stats.norm.ppf(u)\n", " \n", " def Gamma(self, t, xp, u):\n", " return stats.norm.ppf(u, loc=xp)\n", " \n", " def logG(self, t, xp, x):\n", " return np.where((x < self.b) & (x > self.a), 0., -np.inf) \n", " \n", "fk_gp = GaussianProb(a=0., b=1., T=30)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have added: \n", "\n", "* methods `Gamma0` and `Gamma`, which define the deterministic functions $\\Gamma_0$ and $\\Gamma$ we mentioned above. Mathematically, for $U\\sim \\mathcal{U}([0,1]^{d_u})$, then $\\Gamma_0(U)$ is distributed according to $M_0(dx_0)$, and $\\Gamma_t(x_{t-1}, U)$ is distributed according to $M_t(x_{t-1}, dx_t)$. \n", "* class attribute `du`, i.e. $d_u$, the dimension of the $u$-argument of functions $\\Gamma_0$ and $\\Gamma_t$. \n", "\n", "We are now able to run both the SMC and the SQMC algorithms that corresponds to the Feyman-Kac model of interest; let's compare their respective performance. (Recall that function ``multiSMC`` runs several algorithms multiple times, possibly with varying parameters; here we vary parameter `qmc`, which determines whether we run SMC or SMQC.)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAZiUlEQVR4nO3df1AU9/3H8efesbE5o9STINJg0wFMZvJD4sAkEJM0maN/NSZjMwyXH1VGBx2SP1qTqjOKTofYEESd/CxJZEwrFtpxJLadhjo3NDodOh3ONE1MmjRo8xWGSwicEI3mcnj3/cPxmgQsyAKHfF6Pf9jb28/u++4WXnz2c7trxePxOCIiYixXsgsQEZHkUhCIiBhOQSAiYjgFgYiI4RQEIiKGUxCIiBguJdkFjFV3d3eySxARuaxkZmYOO189AhERwykIREQMpyAQETGcgkBExHAKAhERwykIREQMpyAQETHcZXseweWusbGRzs7OZJfBJ598AsC8efOSWkdWVhZ+vz+pNYiYSkFguEgkkuwSRCTJrMv1xjQ6s3h81NTUALBu3bokVyIiE01nFouIyLAUBCIihlMQiIgYTkEgImI4R98a2rNnD0eOHCElJYV58+ZRUVHBzJkzGRwc5OWXX+bYsWO4XC5WrFjBDTfcMOw6Xn/9dVpaWnC73SxevJiHH37YSUkiInKJHAXBzTffzIMPPojb7aahoYHm5mYefvhhAoEAANu3b2dgYIBf/OIXPPXUU7hcX++AHD16lGAwSG1tLbZtMzAw4KQcEREZA0eHhhYtWoTb7QZg4cKFhMNhALq6urjxxhsBSE1NZebMmRw/fnxI+4MHD3Lfffdh23ZiWRERmVzjdkJZa2srRUVFAFx77bUEg0Fuv/12+vr6OH78OL29veTk5HytTSgU4v3336epqQnbtnnkkUeGLHNBIBBI9DSqq6tJS0sbr9KNdiGE9X6KmGvEIKiqqqK/v3/I/NLSUgoKCgDYv38/brebO+64A4C7776brq4uNmzYwNVXX811112X6Dl8VSwW4/Tp02zdupVjx46xc+dOnn/+eSzLGrKsz+fD5/MlHvf29o7+VcpFRaNRQO+niAkudkLZiEFQWVn5P59/4403OHLkCJs3b078AXe73axYsSKxzKZNm5g/f/6Qtl6vl1tvvRXLssjJycHlcnHq1Clmz549UlkiIjJOHB0aeuuttzhw4AA///nPmTFjRmJ+JBIhHo/zrW99i7fffhu3280111wzpH1BQQFHjx7lhhtuoLu7m8HBQWbNmuWkJBEZB1PhoohT5YKIMP0viugoCOrr6xkcHKSqqgqA3NxcysvLGRgYYOvWrbhcLrxeL4899liiTV1dHcXFxWRnZ3PPPffw4osv8vjjj5OSksKjjz467GEhETGPLog4eXTROcPponMyVWnfHH+66JyIiAxLQSAiYjgFgYiI4RQEIiKGUxCIiBhOQSAiYjgFgYiI4RQEIiKGUxCIiBhOQSAiYjgFgYiI4RQEIiKGUxCIiBhOQSAiYjgFgYiI4cbt5vUi4txUuDPYVHHixAngv/clMN1E3iVNQSAyhXR2dvJ+x3/gKm+yS0m+c+fvVvj+xwNJLmQKOB2e0NUrCESmmqu8xBb9MNlVyBTi+ucfJ3b9E7p2ERGZ8hQEIiKGc3RoaM+ePRw5coSUlBTmzZtHRUUFM2fOZHBwkJdffpljx47hcrlYsWIFN9xww5D2H330Ea+88gpffvklbrebVatWkZOT46QkERG5RI56BDfffDPbt2+ntraW+fPn09zcDEAgEABg+/btbNq0iV//+tfEYrEh7RsaGnjggQfYtm0bJSUlNDQ0OClHRETGwFEQLFq0CLfbDcDChQsJh8+PbHd1dXHjjTcCkJqaysyZMzl+/PiQ9pZlcfbsWQDOnDnDnDlznJQjIiJjMG7fGmptbaWoqAiAa6+9lmAwyO23305fXx/Hjx+nt7d3yGGf5cuXs3XrVvbs2UMsFuPJJ5+86PoDgUCip1FdXU1aWtp4lW4027YB9H5OERc+D5Fvsm17wn5PRwyCqqoq+vv7h8wvLS2loKAAgP379+N2u7njjjsAuPvuu+nq6mLDhg1cffXVXHfddYmew1cdPHiQ5cuXc9ttt9HW1kZdXR2VlZXD1uHz+fD5fInHvb29o3uF8j9Fo1FA7+dUceHzEPmmaDTq+Pc0MzNz2PkjBsHF/jBf8MYbb3DkyBE2b96MZZ0/AcTtdrNixYrEMps2bWL+/PlD2h46dIiysjIACgsLeemll0YqR0RExpmjMYK33nqLAwcOsH79embMmJGYH4lE+OKLLwB4++23cbvdXHPNNUPae71e3nvvPQCOHj1KRkaGk3JERGQMHI0R1NfXMzg4SFVVFQC5ubmUl5czMDDA1q1bcblceL1eHnvssUSburo6iouLyc7OZvXq1ezevZtYLIZt26xevdrZqxERkUvmKAiee+65Yeenp6fzzDPPDPvcmjVrEtPXX389Tz/9tJMSRETEIZ1ZLCJiOAWBiIjhFAQiIoZTEIiIGE73IxCZQnp6euD06Qm//rxcZk730dMTmbDVq0cgImI49QhEppD09HTCsRm6Q5l8jeuffyQ9PXXi1j9haxYRkcuCgkBExHAKAhERwykIREQMpyAQETGcgkBExHAKAhERwykIREQMpyAQETGcgkBExHAKAhERwykIREQMpyAQETGco6uPNjU1EQwGsSyL1NRUKioq8Hq9ADQ3N9Pa2orL5aKsrIy8vLwh7U+fPs3OnTv59NNPufrqq/npT3/KVVdd5aQkERG5RI56BEuXLqW2tpZt27axePFi9u3bB0BXVxdtbW3s2LGDjRs3Ul9fTywWG9L+tdde46abbuLZZ5/lpptu4rXXXnNSjoiIjIGjIPB4PInpSCSCZVkAtLe3U1RUhG3bpKenk5GRQUdHx5D27e3t3HXXXQDcddddtLe3OylHRETGwPGNaRobGzl8+DAej4ctW7YAEA6Hyc3NTSzj9XoJh8ND2g4MDDBnzhwA5syZw2effXbR7QQCAQKBAADV1dWkpaU5LV0A27YB9H5OERc+D5Fvsm17wn5PRwyCqqoq+vv7h8wvLS2loKAAv9+P3++nubmZlpYWSkpKiMfj416oz+fD5/MlHvf29o77NkwUjUYBvZ9TxYXPQ+SbotGo49/TzMzMYeePGASVlZWj2sCSJUuorq6mpKSEuXPn0tfXl3guHA4nBpG/KjU1lZMnTzJnzhxOnjzJ7NmzR7UtpxobG+ns7JyUbU11J06cAKCmpibJlUwNWVlZ+P3+ZJchMqkcHRoKhULMnz8fgGAwmEib/Px8nn32WX74wx9y8uRJQqEQOTk5Q9rn5+dz6NAh7r//fg4dOkRBQYGTckats7OT//v3B2S4J2VzU1rKufM/I8c+SG4hU8DH55JdgUhyOAqCvXv3EgqFsCyLtLQ0ysvLgfP/VRUWFrJ27VpcLhcrV67E5To/Ll1XV0dxcTHZ2dncf//97Ny5k9bWVtLS0li7dq3zVzRKGW5YNduatO3J1Lfrs/E/pClyOXAUBE888cRFn1u2bBnLli0bMn/NmjWJ6VmzZrF582YnJYiIiEM6s1hExHAKAhERwykIREQMpyAQETGcgkBExHAKAhERwykIREQMpyAQETGcgkBExHAKAhERwykIREQMpyAQETGcgkBExHAKAhERwykIREQMpyAQETGcgkBExHCO7lAmIhPgdBjXP/+Y7CqS7+xn539eOTu5dUwFp8NA6oStXkEgMoVkZWUlu4Qp48SJAQAWZEzcH8DLR+qE7hsKApEpxO/3J7uEKaOmpgaAdevWJbmS6c9REDQ1NREMBrEsi9TUVCoqKvB6vQA0NzfT2tqKy+WirKyMvLy8Ie337NnDkSNHSElJYd68eVRUVDBz5kwnJYmIyCVyNFi8dOlSamtr2bZtG4sXL2bfvn0AdHV10dbWxo4dO9i4cSP19fXEYrEh7W+++Wa2b99ObW0t8+fPp7m52Uk5IiIyBo6CwOPxJKYjkQiWZQHQ3t5OUVERtm2Tnp5ORkYGHR0dQ9ovWrQIt9sNwMKFCwmHw07KERGRMXA8RtDY2Mjhw4fxeDxs2bIFgHA4TG5ubmIZr9c74h/51tZWioqKLvp8IBAgEAgAUF1dTVpa2phrtm2byJhby3Rm27ajfUvGj23bAPo8JsGIQVBVVUV/f/+Q+aWlpRQUFOD3+/H7/TQ3N9PS0kJJSQnxePySiti/fz9ut5s77rjjosv4fD58Pl/icW9v7yVt46ui0eiY28r0Fo1GHe1bMn4u/J7q8xg/mZmZw84fMQgqKytHtYElS5ZQXV1NSUkJc+fOpa+vL/FcOBxODCJ/0xtvvMGRI0fYvHlz4tCSiIhMHkdjBKFQKDEdDAYTaZOfn09bWxvRaJSenh5CoRA5OTlD2r/11lscOHCA9evXM2PGDCeliIjIGDkaI9i7dy+hUAjLskhLS6O8vBw4f1JMYWEha9euxeVysXLlSlyu85lTV1dHcXEx2dnZ1NfXMzg4SFVVFQC5ubmJdYiIyORwFARPPPHERZ9btmwZy5YtGzJ/zZo1iennnnvOyeZFRGQc6KJzIiKGUxCIiBhOQSAiYjgFgYiI4RQEIiKGUxCIiBhOQSAiYjgFgYiI4RQEIiKGUxCIiBhOQSAiYjgFgYiI4RQEIiKGUxCIiBhOQSAiYjgFgYiI4RzdmEZEpqfGxkY6OzuTWsOJEycAqKmpSWodcP6ui36/P9llTBgFgYhMSbqP+eRREIjIENP5v18ZSmMEIiKGc9QjaGpqIhgMYlkWqampVFRU4PV6AWhubqa1tRWXy0VZWRl5eXkXXc/vf/97Ghoa2LVrF7Nnz3ZSkoiIXCJHPYKlS5dSW1vLtm3bWLx4Mfv27QOgq6uLtrY2duzYwcaNG6mvrycWiw27jt7eXt555x3S0tKclCIiImPkKAg8Hk9iOhKJYFkWAO3t7RQVFWHbNunp6WRkZNDR0THsOn71q1/x0EMPJdqKiMjkcjxY3NjYyOHDh/F4PGzZsgWAcDhMbm5uYhmv10s4HB7SNhgM4vV6ufbaa0fcTiAQIBAIAFBdXe2oB2HbNpExt5bpzLZt9U7FOCMGQVVVFf39/UPml5aWUlBQgN/vx+/309zcTEtLCyUlJcTj8RE3HIlE2L9/P5s2bRpVoT6fD5/Pl3jc29s7qnbDiUajY24r01s0GnW0b4lMZZmZmcPOHzEIKisrR7WBJUuWUF1dTUlJCXPnzqWvry/xXDgcTgwiX/DJJ5/Q09PDz372MwD6+vpYv349Tz31FN/+9rdHtU0REXHO0RhBKBRKTAeDwUTa5Ofn09bWRjQapaenh1AoRE5OztfaLliwgF27dvHCCy/wwgsvMHfuXJ5++mmFgIjIJHM0RrB3715CoRCWZZGWlkZ5eTlw/nTswsJC1q5di8vlYuXKlbhc5zOnrq6O4uJisrOznVcvIiKOWfHRHNCfgrq7u8fctqamhsixD1g1W99Ukv/a9VmcGdnXsW7dumSXIjIhLjZGoDOLRUQMpyAQETGcgkBExHAKAhERwykIREQMZ+T9CHp6ejgzeP5bIiIXhAbB09OT7DJEJp16BCIihjOyR5Cenk7k1EmdRyBfs+uzODPS05NdhsikU49ARMRwCgIREcMpCEREDKcgEBExnIJARMRwCgIREcMpCEREDKcgEBExnIJARMRwCgIREcMpCEREDKcgEBExnKOLzjU1NREMBrEsi9TUVCoqKvB6vQA0NzfT2tqKy+WirKyMvLy8Ydfx+uuv09LSgtvtZvHixTz88MNOShIRkUvkKAiWLl1KaWkpAH/605/Yt28f5eXldHV10dbWxo4dOzh58iRVVVU888wzuFxf74AcPXqUYDBIbW0ttm0zMDDgpBwRERkDR4eGPB5PYjoSiWBZ5y/r3N7eTlFREbZtk56eTkZGBh0dHUPaHzx4kPvuuw/btgFITU11Uo6IiIyB4/sRNDY2cvjwYTweD1u2bAEgHA6Tm5ubWMbr9RIOh4e0DYVCvP/++zQ1NWHbNo888gg5OTnDbicQCBAIBACorq4mLS1tzDXbtk1kzK1lOrNt29G+JXI5GjEIqqqq6O/vHzK/tLSUgoIC/H4/fr+f5uZmWlpaKCkpIR4f3S0gY7EYp0+fZuvWrRw7doydO3fy/PPPJ3oWX+Xz+fD5fInHvb29o9rGcKLR6JjbyvQWjUYd7VsiU1lmZuaw80cMgsrKylFtYMmSJVRXV1NSUsLcuXPp6+tLPBcOhxODyF/l9Xq59dZbsSyLnJwcXC4Xp06dYvbs2aPapoiIOOdojCAUCiWmg8FgIm3y8/Npa2sjGo3S09NDKBQa9pBPQUEBR48eBaC7u5vBwUFmzZrlpCQREblEjsYI9u7dSygUwrIs0tLSKC8vByArK4vCwkLWrl2Ly+Vi5cqViW8M1dXVUVxcTHZ2Nvfccw8vvvgijz/+OCkpKTz66KPDHhYSEZGJY8VHe0B/iunu7h5z25qaGv7v3x+Q4R7Hgi5TfefO/5yr94KPz8F3F17HunXrkl2KyIQY8xjBdJSVlZXsEqaMwRMnAJixYEGSK0m+76J9Q8xkZI9A/qumpgZA/wWLGOBiPQJda0hExHAKAhERwykIREQMpyAQETGcgkBExHAKAhERwykIREQMpyAQETGcgkBExHAKAhERwykIREQMpyAQETGcgkBExHAKAhERwykIREQMpyAQETGcgkBExHAKAhERwzm6Z3FTUxPBYBDLskhNTaWiogKv1wtAc3Mzra2tuFwuysrKyMvLG9L+o48+4pVXXuHLL7/E7XazatUqcnJynJQkIiKXyFGPYOnSpdTW1rJt2zYWL17Mvn37AOjq6qKtrY0dO3awceNG6uvricViQ9o3NDTwwAMPsG3bNkpKSmhoaHBSjoiIjIGjIPB4PInpSCSCZVkAtLe3U1RUhG3bpKenk5GRQUdHx5D2lmVx9uxZAM6cOcOcOXOclCMiImPg6NAQQGNjI4cPH8bj8bBlyxYAwuEwubm5iWW8Xi/hcHhI2+XLl7N161b27NlDLBbjySefvOh2AoEAgUAAgOrqatLS0pyWLoBt2wB6P0UMNmIQVFVV0d/fP2R+aWkpBQUF+P1+/H4/zc3NtLS0UFJSQjweH9XGDx48yPLly7nttttoa2ujrq6OysrKYZf1+Xz4fL7E497e3lFtQ/63aDQK6P0UMUFmZuaw80cMgov9Yf6mJUuWUF1dTUlJCXPnzqWvry/xXDgcTgwif9WhQ4coKysDoLCwkJdeemlU2xIRkfHjaIwgFAolpoPBYCJt8vPzaWtrIxqN0tPTQygUGvbbQF6vl/feew+Ao0ePkpGR4aQcEREZA0djBHv37iUUCmFZFmlpaZSXlwOQlZVFYWEha9euxeVysXLlSlyu85lTV1dHcXEx2dnZrF69mt27dxOLxbBtm9WrVzt/RSIickms+GgP6E8x3d3dyS5hWqipqQFg3bp1Sa5ERCbaxcYIdGaxiIjhFAQiIoZTEIiIGE5BICJiOAWBiIjhFAQiIoZTEIiIGE5BICJiOJ1QliSNjY10dnYmuwxOnDgBwIIFC5JaR1ZWFn6/P6k1iEx3Y77onExvM2bMSHYJIpJk6hGIiBhCl5gQEZFhKQhERAynIBARMZyCQETEcAoCERHDKQhERAynIBARMZyCQETEcJftCWUiIjI+1CMQNmzYkOwSRIalfXNyKAhERAynIBARMZyCQPD5fMkuQWRY2jcnhwaLRUQMpx6BiIjhFAQiIoZTEIiIGE5BICJiON2zeBr74osv2LlzJ+FwmFgsxo9+9CP27t3L7bffzrvvvsu5c+coLy+nsbGRjz/+mHvvvZcf/OAHABw4cIDDhw/jcrnIy8vjoYceSvKrkcvdcPujx+Ph1VdfZdasWXzve9+jp6eHDRs28Lvf/Y6enh76+/sJhUL8+Mc/5sMPP+Qf//gHXq+X9evXk5KSQkdHB6+++iqRSISUlBQ2b97MlVdemeyXevmJy7T1t7/9Lf7LX/4y8fjzzz+PV1RUxP/85z/H4/F4fPfu3fHHH388fubMmfjAwEB85cqV8Xg8Hn/zzTfjGzdujH/xxRfxeDweP3Xq1OQXL9POcPvjmjVr4t3d3fFYLBbfvn17/KmnnorH4/H4b3/72/imTZvi0Wg0/p///Cf+0EMPxd988814PB6P19TUxP/+97/Ho9Fo/NFHH41/+OGHifUNDg5O/gubBnRoaBpbsGAB77zzDg0NDfzrX//C4/EAkJ+fn3g+JyeHK6+8ktmzZ2PbNp9//jnvvPMO3//+95kxYwYAV111VdJeg0wf39wfe3p6SE9PZ/78+ViWxZ133vm15W+55RZSUlJYsGABsViMvLy8xHo+/fRTuru7mTNnDjk5OQB4PB7cbvekv67pQIeGprHMzEyefvpp3nzzTX7zm9+waNEiAFJSzn/sLpcL27YTy7tcLs6dO0c8HseyrKTULNPXxfbHi/nqfup2uxP7pGVZif1Uxod6BNNYOBzmiiuu4M477+Tee+/l+PHjo2q3aNEi/vKXvxCJRAA4ffr0RJYphvjm/vjBBx/Q09PDxx9/DMBf//rXS1rfd77zHU6ePElHRwcAZ8+e5dy5c+NetwnUI5jGTpw4QUNDA5ZlkZKSwqpVq9ixY8eI7fLy8vjoo4/YsGEDKSkp3HLLLTz44IOTULFMZ8Ptj6dOnaK6uppZs2Zx/fXX09nZOer1paSk8JOf/ITdu3fz5ZdfcsUVV1BZWanDQ2OgS0yIyJTw7rvv8oc//EGXnk4CHRoSETGcegQiIoZTj0BExHAKAhERwykIREQMpyAQETGcgkBExHD/D4az9ol+emy4AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "results = particles.multiSMC(fk=fk_gp, qmc={'smc': False, 'sqmc': True}, N=100, nruns=10)\n", "\n", "sb.boxplot(x=[r['qmc'] for r in results], y=[r['output'].logLt for r in results]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We do get some variance reduction, but not so much. Let's see if we can do better by addressing point 1 above.\n", "\n", "The considered problem has the structure of a state-space model, where process $(X_t)$ is a random walk, $Y_t = \\mathbb{1}_{[a,b]}(X_t)$, and $y_t=1$ for all $t$'s. This remark leads us to define alternative Feynman-Kac models, that would correspond to *guided* and *auxiliary* formalisms of that state-space model. In particular, for the guided filter, the optimal proposal distribution, i.e. the distribution of $X_t|X_{t-1}, Y_t$, is simply a Gaussian distribution truncated to interval $[a, b]$; let's implement the corresponding Feynman-Kac class. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def logprobint(a, b, x):\n", " \"\"\" returns log probability that X_t\\in[a,b] conditional on X_{t-1}=x\n", " \"\"\"\n", " return np.log(stats.norm.cdf(b - x) - stats.norm.cdf(a - x))\n", "\n", "class Guided_GP(GaussianProb):\n", "\n", " def Gamma(self, t, xp, u):\n", " au = stats.norm.cdf(self.a - xp)\n", " bu = stats.norm.cdf(self.b - xp)\n", " return xp + stats.norm.ppf(au + u * (bu - au))\n", "\n", " def Gamma0(self, u):\n", " return self.Gamma(0, 0., u)\n", "\n", " def M(self, t, xp):\n", " return self.Gamma(t, xp, stats.uniform.rvs(size=xp.shape))\n", "\n", " def M0(self, N):\n", " return self.Gamma0(stats.uniform.rvs(size=N))\n", "\n", " def logG(self, t, xp, x):\n", " if t == 0:\n", " return np.full(x.shape, logprobint(self.a, self.b, 0.))\n", " else:\n", " return logprobint(self.a, self.b, xp)\n", " \n", "fk_guided = Guided_GP(a=0., b=1., T=30) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this particular case, it is a bit more convenient to define methods `Gamma0` and `Gamma` first, and then define methods `M0` and `M`. \n", "\n", "To derive the APF version, we must define the auxiliary functions (functions $\\eta_t$ in Chapter 10 of the book) that modify the resampling probabilities; in practice, we define the log of these functions, \n", "as follows:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "class APF_GP(Guided_GP): \n", " def logeta(self, t, x):\n", " return logprobint(self.a, self.b, x)\n", " \n", "fk_apf = APF_GP(a=0., b=1., T=30) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok, now everything is set! We can do a 3x2 comparison of SMC versus SQMC, for the 3 considered Feynman-Kac models. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3df1RUdf4/8OfMMJBADOAAamIqYvaD8pi0aRh6gtqTwbHWH+DRsuPGuqN1VrSwk7/2w1a0/uBYW5FldcTfuoG5uXoiA0ItfrSe0mwT1EQdRZxhAJFxZnh//3C9X234MXBhBuY+H/84c+e+7329544859733HtVQggBIiJSLLWnCyAiIs9iEBARKRyDgIhI4RgEREQKxyAgIlI4BgERkcL5eLqArjp//rynSyAi6lMGDRrU6nTuERARKRyDgIhI4RgEREQKxyAgIlI4BgERkcIxCIiIFI5BQESkcH32PAI5tmzZgurqapfmvXjxIgAgIiLCpfkjIyMxc+bMLtdGRORuigyCzrBarZ4ugYioR6n66o1p3HVm8VtvvQUAyMjIcMv6iIh6Cs8sJiKiVjEIiIgUjkFARKRwDAIiIoVjEBARKRyDgIhI4RgEREQKxyAgIlI4WWcW5+bmoqKiAj4+PoiIiIDBYEBAQADsdjtycnJw6tQptLS04NFHH8XTTz/t1H7Hjh346quvEBQUBABITU3FmDFj5JRERESdJCsI7r//fsycORMajQabNm1CXl4eZs2ahW+//RZ2ux1r1qyB1WpFeno6HnnkEYSHhzstY/LkyUhOTpZTBhERySDr0NADDzwAjUYDABg5ciRMJpP0WnNzMxwOB65duwYfHx/4+/vLq5SIiHpEt1107sCBAxg/fjwA4OGHH0Z5eTnS0tJw7do1PPfccwgMDGy13f79+1FcXIzhw4fj2WefbXO+goICFBQUAACysrKg1+u7q/R2abVaAHDb+oiI3K3DIMjMzERdXZ3T9JSUFMTGxgIAPvvsM2g0GkyYMAEAUFlZCbVajQ8++ABXrlzB8uXLERMT43Qp58cffxxTp04FAGzfvh0bN26EwWBotY6EhAQkJCRIz2tra13sojw2m82t6yMi6iltXXSuwyBYtmxZu68XFhaioqICy5cvh0qlAgCUlJRg9OjR8PHxgU6nw1133YWqqiqnIAgODpYeP/bYY9KVPomIyH1kjREcOXIEu3fvRkZGBvz8/KTper0eR48ehRACzc3NOHHiBO644w6n9mazWXpcWlqKyMhIOeUQEVEXyBoj2LBhA+x2OzIzMwEA0dHRSEtLw+9//3u89957WLRoEYQQmDRpEu68804AQE5ODhITExEVFYVNmzbh9OnTUKlUCAsLQ1pamvweERFRp/DGNB3gjWmIyFvwxjRERNQqBgERkcIxCIiIFI5BQESkcAwCIiKFYxAQESkcg4CISOEYBERECscgICJSOAYBEZHCMQiIiBSOQUBEpHAMAiIihWMQEBEpHIOAiEjhGARERArHICAiUjhZt6rMzc1FRUUFfHx8EBERAYPBgICAANjtdqxfvx5VVVVQq9WYM2cO7r33Xqf2jY2NyM7OxqVLlxAWFoaFCxciMDBQTklERNRJsvYI7r//fqxZswarV6/GwIEDkZeXBwAoKCgAAKxZswZLly7Fxo0b0dLS4tQ+Pz8fMTExePvttxETE4P8/Hw55RARURfICoIHHngAGo0GADBy5EiYTCYAwNmzZ3HfffcBAHQ6HQICAnDy5Emn9mVlZYiPjwcAxMfHo6ysTE45RETUBbIODd3swIEDGD9+PABg6NChKC8vxyOPPILLly/j5MmTqK2txYgRI25pY7FYEBISAgAICQlBfX19m8svKCiQ9jSysrKg1+u7q/R2abVaAHDb+tzNZDJh1apVeOWVV6RtQUTK0mEQZGZmoq6uzml6SkoKYmNjAQCfffYZNBoNJkyYAACYNGkSzp49iyVLliAsLAx33XWXtOfQVQkJCUhISJCe19bWylqeq2w2m1vX5265ubn46aef8Omnn2L27NmeLoeIetCgQYNand5hECxbtqzd1wsLC1FRUYHly5dDpVIBADQaDebMmSPNs3TpUgwcONCprU6ng9lsRkhICMxmM4KCgjoqh7pRXV0dSkpKIIRASUkJkpOTodPpPF0WEbmZrDGCI0eOYPfu3cjIyICfn5803Wq1orm5GQDwww8/QKPRYPDgwU7tx44di6KiIgBAUVGRtIdB7rFnzx5pEL+lpQWff/65hysiIk+QNUawYcMG2O12ZGZmAgCio6ORlpYGi8WC119/HWq1GqGhoViwYIHUJicnB4mJiYiKisKUKVOQnZ2NAwcOQK/XIz09XV5vqFMOHz4Mh8MBAHA4HDh8+DAPDxEpkEoIITxdRFecP3/eLet56623AAAZGRluWZ875ebmori4GA6HAxqNBo8++iiDgMiLtTVGwDOLFSwpKQlq9fWPgFqtRnJysocrIiJPYBAoWHBwMOLi4qBSqRAXF8eBYiKF6rbzCKhvSkpKwrlz57g3QKRgDAKFCw4OxpIlSzxdBhF5EA8NEREpHIOAiEjhGARERArHICAiUjgGARGRwjEIiIgUjkFARKRwDAIiIoVjEBARKRyDgIhI4RgEREQKx2sNeZktW7agurra5fkvXrwIAIiIiHBp/sjISMycObNLtRFR7yQrCHJzc1FRUQEfHx9ERETAYDAgICAAdrsd69evR1VVFdRqNebMmYN7773Xqf2OHTvw1VdfSfcqTk1NxZgxY+SURJ1ktVo9XQIReZisILj//vsxc+ZMaDQabNq0CXl5eZg1axYKCgoAAGvWrIHFYsEbb7yBN998U7oJys0mT57MSyB3o85+W/fmO7ARkWtkjRE88MAD0Gg0AICRI0fCZDIBAM6ePYv77rsPAKDT6RAQEICTJ0/KLJWIiHpCtw0WHzhwAKNHjwYADB06FOXl5XA4HKipqcHJkydRW1vbarv9+/dj8eLFeO+999DY2Nhd5RARkYs6PDSUmZmJuro6p+kpKSmIjY0FAHz22WfQaDSYMGECAGDSpEk4e/YslixZgrCwMNx1113SnsPNHn/8cUydOhUAsH37dmzcuBEGg6HVOgoKCqRDTllZWdDr9S52UR6tVgsAblufu3l7/4ioYx0GwbJly9p9vbCwEBUVFVi+fDlUKhUAQKPRYM6cOdI8S5cuxcCBA53aBgcHS48fe+wx6Xh1axISEpCQkCA9b2sPo7vZbDa3rs/dvL1/RPT/DRo0qNXpsg4NHTlyBLt370ZGRgb8/Pyk6VarFc3NzQCAH374ARqNBoMHD3ZqbzabpcelpaWIjIyUUw4REXWBrF8NbdiwAXa7HZmZmQCA6OhopKWlwWKx4PXXX4darUZoaCgWLFggtcnJyUFiYiKioqKwadMmnD59GiqVCmFhYUhLS5PXGyIi6jRZQfDOO++0Oj08PBzr1q1r9bV58+ZJj1988UU5qyciom7AS0wQESkcg4CISOEYBERECscgICJSOK+5+mhnr7rpqjNnzgBAu+c4dBWv5ElEvYHXBEF1dTV+/eW/GOh8ArMsWsf1f69V/bdbl2t0dOviiIi6zGuCAAAGaoA/6vrG0a6PLC2eLoGICADHCIiIFI9BQESkcAwCIiKF86oxAm/VU7+IAvirKCJiEPQJ1dXV+LnyFBDYv/sX7rh+6fCfL9R373IbL3fv8oioxzAI+orA/hCjn/J0FS5THfmXp0tAXV0dcnJy8Oc//xk6nc7T5VAncfu5D8cIyGvt2rULv/zyC3bu3OnpUqgL9uzZgxMnTuDzzz/3dClej0FAXqmurg6HDx8GABw+fBgWi8XDFXW/uro6ZGVleW3fSkpKIIRASUmJV/bx2LFjmDt3Ln766SdPl8IgIO+0a9cuCCEAAEIIr9wr8OZvzHv27EFLy/WTLltaWryyj++//z6EEHjvvfc8XQrHCKhvcfUXVP/9762XBDl06BAuX+54ALuv/NKprq4O33zzjfSNOTk5uU8cR3d1+504cUIKAofDgaKiIpw/f77dNn1l2wHX9waampoAAE1NTfjpp59wzz33eKweWUGwbds2lJeXQ6VSQafTwWAwIDQ0FACQl5eHAwcOQK1W4/nnn8fo0aOd2jc2NiI7OxuXLl1CWFgYFi5ciMDAQDklUR/UmZ/H1tTUSPfD7qwbP5XtaPmu1tITf3hcfS8uXLgAu90OALDZbFi5ciUGDBjQbpue+kPZE9tPpVI5Pe9o+3l62wGdC7qbrV27FtHR0e226cmgkxUEycnJSElJAQDs3bsXu3btQlpaGs6ePYtDhw5h7dq1MJvNyMzMxLp166BW33okKj8/HzExMZgyZQry8/ORn5+PWbNmySmJ+qCKigrUmevg4+PXbctUqdQQouWW53ab6LBdo+0qGht+7XA+u92Kmpoal/5jbtmyBQcPHuxwPuD6H/Ub34Tb89t5LBYLGhoa2m1TVVXlch0A8Mgjj7jUv4qKCljqzPDz7b4jzVoN4HDc/FwALdZ22zQ3WXHmdMdjCdZrLS5vO8A926+lpcUpHH6rM9vP1W13g6wt5+/vLz22Wq1SipeVlWH8+PHQarUIDw/HgAEDUFlZ6dS+rKwM8fHxAID4+HiUlZXJKYdIotFo231OvdvNOwSq3zyn7id7jGDr1q0oLi6Gv78/VqxYAQAwmUy37OaEhobCZDI5tbVYLAgJCQEAhISEoL6+7ZOaCgoKUFBQAADIysqCXq+/5XWtVotrcjvjZlqt1qkfbc3XF7nav7i4OJw6dcqlZRqNRly9etWlef935AQajQa+vq5dn7xfv34YOHCgS/MOGzbMpf699NJLeOmll1xa5ocffujSe3Hs2DFpMBy4fuiko2PMw4YNwwsvvOBSHZ3RE9vvet+uzycAQO3ndLjot3pi2wE9s/2OHj3qNM1T2w9wIQgyMzNRV1fnND0lJQWxsbFITU1Famoq8vLysG/fPkyfPv2WD2h3SUhIQEJCgvS8trb2ltdtNlu3r7On2Ww2p3605ty5c0DjlV5xkpbLGi/j3Llml/r39NNPu7zYzhyPNhqNuHLlCoYOHQofH9e+83T2OKwr/esMV9+L3NxcfP3119LziRMnYvbs2R226+56gZ7ZfhcuXLglMHx9fbt9DMST78Vvt9+kSZPcsv0GDRrU6vQO/3csW7bMpRXExcUhKysL06dPR//+/W/5hYbJZJIGkW+m0+lgNpsREhICs9mMoKAgl9ZFytVXfhXS05KSkvDNN9/AbrfDx8cHycnJni7JJa5uP4PBcMtzq9WKjIyMnijJI5KSklBUVISWlhao1WqPbz9Zh4aMRqO0K1ZeXi6lzdixY/H222/jqaeegtlshtFoxIgRI5zajx07FkVFRZgyZQqKiooQGxsrpxyvFR4eDlNLfZ+7xER4OIO9pwQHB2PChAkoLCzEhAkT+sRPRztj3LhxKC4uhsPhgEajwbhx4zxdUrcKDg5GfHw8CgsLER8f7/HtJysINm/eDKPRCJVKBb1ej7S0NADXd9HGjRuH9PR0qNVqzJ07V/rFUE5ODhITExEVFYUpU6YgOzsbBw4cgF6vR3p6uvweESlEUlISzp075/Fvkz0hKSkJJSUlcDgcveIbc0/oTdtPVhAsXry4zdeeeeYZPPPMM07T582bJz2+/fbbsXz5cjklEClWcHAwlixZ4ukyekRwcDDi4uJQWFiIuLg4j39j7gm9afvxzGIi6pV60zdmb8cgIKJeqTd9Y/Z2vOgcEZHCMQiIiBSOh4b6isbLPXNC2dX/XZulXzcPxjVeBsCfjxL1BQyCPiAyMrLHln3mzPWzxocM6O4/2kE9WjcRdR8GQR/Qk2fTvvXWWwDgVWdtElHncIyAiEjhvGaPoKamBlftwEeWjq8F3hsY7UC/mhpPl0FExD0CIiKl85o9gvDwcFxrMOOPur6RbR9ZWuAbHu7pMoiIuEdARKR0DAIiIoVjEBARKRyDgIhI4RgEREQKxyAgIlI4WT8f3bZtG8rLy6FSqaDT6WAwGKSb1Ofl5eHAgQNQq9V4/vnnMXr0aKf2O3bswFdffSXdtD41NRVjxoyRUxIREXWSrCBITk5GSkoKAGDv3r3YtWsX0tLScPbsWRw6dAhr166F2WxGZmYm1q1bJ923+GaTJ0/mHYiIiDxI1qEhf39/6bHVaoVKpQIAlJWVYfz48dBqtQgPD8eAAQNQWVkpr1IiIuoRss8s3rp1K4qLi+Hv748VK1YAAEwmE6Kjo6V5QkNDYTKZWm2/f/9+FBcXY/jw4Xj22WcRGBjY6nwFBQUoKCgAAGRlZUGv19/yularxTW5nXEzrVbr1A9P1ADA43UQked0GASZmZmoq6tzmp6SkoLY2FikpqYiNTUVeXl52LdvH6ZPnw4hhEsrf/zxxzF16lQAwPbt27Fx40YYDIZW501ISEBCQoL0vLa29pbXbTabS+vsTWw2m1M/PFED4Px+EpH3GTRoUKvTOwyCZcuWubSCuLg4ZGVlYfr06ejfvz8uX74svWYymaRB5JsFBwdLjx977DHp2vhEROQ+ssYIjEaj9Li8vFxKm7Fjx+LQoUOw2WyoqamB0WjEiBEjnNqbzWbpcWlpKe9oRUTkAbLGCDZv3gyj0QiVSgW9Xo+0tDQA12+tOG7cOKSnp0OtVmPu3LnSL4ZycnKQmJiIqKgobNq0CadPn4ZKpUJYWJjUnoiI3EdWECxevLjN15555hk888wzTtPnzZsnPX7xxRflrJ6IiLoBzywmIlI4BgERkcIxCIiIFI5BQESkcF5zz2K6bsuWLaiurnZ5/jNnzgCAy+dwREZGYubMmV2qjYh6J68KAqPj+k3hu9Nlx/V/+2u6dbEwOoA7u3eRXeLn5+fpEojIw7wmCHrqZDTb/74x+w4Z0q3LvRM9UzO/rRNRZ6mEqxcG6mXOnz/vlvXcOGSSkZHhlvUREfWUtq41xMFiIiKFYxAQESkcg4CISOEYBERECscgICJSOAYBEZHCMQiIiBSOQUBEpHAMAiIihZN1iYlt27ahvLwcKpUKOp0OBoMBoaGhaGhowNq1a1FZWYmJEydi7ty5rbZvbGxEdnY2Ll26hLCwMCxcuBCBgYFySiIiok6StUeQnJyM1atXY9WqVRgzZgx27doFANBqtZgxYwZmz57dbvv8/HzExMTg7bffRkxMDPLz8+WUQ11QV1eHrKwsWCwWT5dCRB4iKwj8/f2lx1arFSqVCgBw2223YdSoUfD19W23fVlZGeLj4wEA8fHxKCsrk1MOdcGePXtw4sQJfP75554uhYg8RPbVR7du3Yri4mL4+/tjxYoVnWprsVgQEhICAAgJCUF9fX2b8xYUFKCgoAAAkJWVBb1e3/WiO0Gr1QKA29bnTiaTCQcPHoQQAgcPHsScOXOk7UFEytFhEGRmZqKurs5pekpKCmJjY5GamorU1FTk5eVh3759mD59eo8UmpCQgISEBOl5bW1tj6znt2w2m1vX5065ublwOK7fcMHhcODTTz/t8HAeEfVdbV19tMMgWLZsmUsriIuLQ1ZWVqeCQKfTwWw2IyQkBGazGUFBQS63JfkOHz58SxAcPnyYQUCkQLLGCIxGo/S4vLy8zbRpy9ixY1FUVAQAKCoqQmxsrJxyqJPGjRsHjeb6rdc0Gg3GjRvn4YqIyBNkjRFs3rwZRqMRKpUKer0eaWlp0mvz589HU1MT7HY7ysrKsHTpUgwePBg5OTlITExEVFQUpkyZguzsbBw4cAB6vR7p6emyO0SuS0pKQklJCRwOB9RqNZKTkz1dEhF5AO9Q1gFvv0NZbm4uCgsLMXHiRB4WIvJyXR4jIO+WlJSEc+fOcW+ASMEYBAoXHByMJUuWeLoMIvIgXmuIiEjhGARERArHICAiUjgGARGRwjEIiIgUjkFARKRwDAIiIoVjEBARKRyDgIhI4RgEREQKxyAgIlI4BgERkcIxCIiIFI5BQESkcLIuQ71t2zaUl5dDpVJBp9PBYDAgNDQUDQ0NWLt2LSorKzFx4kTMnTu31fY7duzAV199Jd2rODU1FWPGjJFTEhERdZKsIEhOTkZKSgoAYO/evdi1axfS0tKg1WoxY8YMnDlzBtXV1e0uY/LkybwpChGRB8k6NOTv7y89tlqtUKlUAIDbbrsNo0aNgq+vr7zqiIiox8m+Q9nWrVtRXFwMf39/rFixotPt9+/fj+LiYgwfPhzPPvssAgMDW52voKAABQUFAICsrCzo9XpZdbtKq9UCgNvWR0Tkbh3evD4zMxN1dXVO01NSUhAbGys9z8vLg81mw/Tp06VphYWFqKqqanOMoK6uThof2L59O8xmMwwGg0uF8+b1RESd0+Wb1y9btsylFcTFxSErK+uWIOhIcHCw9Pixxx6T/ugSEZH7yBojMBqN0uPy8vI206YtZrNZelxaWorIyEg55RARURfIGiPYvHkzjEYjVCoV9Ho90tLSpNfmz5+PpqYm2O12lJWVYenSpRg8eDBycnKQmJiIqKgobNq0CadPn4ZKpUJYWNgt7YmIyD06HCPorThGQETUOW0dteGZxURECscgICJSOAYBEZHCMQiIiBSOQUBEpHAMAiIihWMQEBEpHIOAiEjhGARERArHICAiUjgGARGRwjEIiIgUjkFARKRwDAIiIoVjEBARKRyDgIhI4RgEREQKJ+tWldu2bUN5eTlUKhV0Oh0MBgNCQ0Pxww8/YPPmzbDb7fDx8cHs2bNx3333ObVvbGxEdnY2Ll26hLCwMCxcuBCBgYFySiIiok6SdavKpqYm+Pv7AwD27t2Ls2fPIi0tDadOnYJOp0NoaCjOnDmD119/HR988IFT+02bNiEwMBBTpkxBfn4+GhsbMWvWLJfWzVtVEhF1To/cqvJGCACA1WqFSqUCAAwbNgyhoaEAgMjISNhsNthsNqf2ZWVliI+PBwDEx8ejrKxMTjlERNQFsg4NAcDWrVtRXFwMf39/rFixwun17777DsOGDYNWq3V6zWKxICQkBAAQEhKC+vr6NtdTUFCAgoICAEBWVhb0er3c0l1yo253rY+IyN06DILMzEzU1dU5TU9JSUFsbCxSU1ORmpqKvLw87Nu3D9OnT5fmqa6uxubNm/Haa6/JLjQhIQEJCQnS89raWtnLdMWNPRl3rY+IqKe0dWiowyBYtmyZSyuIi4tDVlaWFASXL1/G6tWrMX/+fAwYMKDVNjqdDmazGSEhITCbzQgKCnJpXURE1H1kjREYjUbpcXl5uZQ2V65cQVZWFlJTUzFq1Kg2248dOxZFRUUAgKKiIsTGxsoph4iIukDWr4ZWr14No9EIlUoFvV6PtLQ0hIaG4p///Cfy8/Nv2RNYunQpdDodcnJykJiYiKioKDQ0NCA7Oxu1tbXQ6/VIT093+eejcn41tGXLFlRXV7s075kzZwAAQ4YMcWn+yMhIzJw5s8u1ERH1lLYODckKAk9yVxBcvHgRABAREeHS/AwCIuqtGARERArXI+cREBFR38cgICJSOAYBEZHCMQiIiBSOQUBEpHAMAiIihWMQEBEpHIOAiEjh+uwJZURE1D24R+CCJUuWeLqEHuXN/fPmvgHsX1/XW/rHICAiUjgGARGRwmlWrly50tNF9AXDhw/3dAk9ypv75819A9i/vq439I+DxURECsdDQ0RECscgICJSOK8KgpqaGixatEj2cr744gtYrdZuqEgeb+tPV1RVVeHjjz9u9bX58+ejvr7e5WUVFhZiw4YN3VVal3hbf7rDuXPn8PLLL+OVV17BhQsXPF1Ot+hrffLxdAG90d69ezFhwgT4+fl5upRu0Zf7ExUVhaioKE+X0W28rT/doaysDLGxsZg+fbqnS+k2fa1PXhcEDocD//jHP3D69GkMHDgQCxYswC+//ILc3Fw4HA5ERUXhhRdegFarxY8//ug0/csvv4TJZMJf//pXBAUFYcWKFbcsv7q6Gu+99x7sdjuEEFi0aBE0Gg3eeOMNjBo1CidOnMCdd96JiRMnYufOnbBYLHjppZcwYsQINDc34+OPP0ZVVRVUKhWmTp2Khx9+uNf1Z+DAgfjss89QVFQEvV6P22+/HcOHD0dycjJWrlyJoUOH4tSpU6ivr8f8+fORn5+PM2fOYPz48UhJSQEAFBUVYc+ePVCpVBgyZAhefPFFaZ27du1CSUkJ+vfvLy37+++/x+zZsxEVFYX6+nq8+uqrePfdd3Hs2DHs2bMHS5YsQUNDA9atW4f6+nqMGDECN//Oobi4GP/+979ht9sRHR2NP/7xj1Cr1fj666+Rn5+P4OBgDBw4EFqt1uk9bmlpwfvvv4+TJ08CACZNmoSnnnrK5b6uXr0a33//PXx8fHD77bfjiSee6JX9OXnyJN5//334+vpi1KhROHLkCNasWYPCwkKUlpaipaUF1dXVSEpKgt1uR3FxMbRaLV599VUEBgbiwoUL+PDDD1FfXw+1Wo2FCxdiwIABrX5u//73v+Py5cuw2Wx48sknkZCQgNmzZyMxMRHHjh1DQEAA/vKXv6CyshJffPEF1Go1jh8/fsvnU+52ae8z2BU92Sd3bptWCS9y8eJFMW3aNHH8+HEhhBDvvvuu2LVrl5g3b544d+6cEEKId955R/zrX/8SVqu11elCCGEwGITFYml1HRs2bBDFxcVCCCFsNpuwWq3i4sWLYsaMGeLXX38VDodDvPLKK+Ldd98VLS0torS0VLz11ltCCCFyc3PFJ598Ii2roaGhV/anqqpKpKeni+bmZnHlyhWxYMECsXv3biGEECtWrBC5ublCCCG++OILkZaWJkwmk7h27Zr405/+JOrr68WZM2fESy+9JK3z5n5WVlaKxYsXC6vVKpqamsSLL74odu/eLVasWCEqKyuFEEJYLBZhMBiEEEIcPXpUvPnmm1KtO3fuFEIIUVFRIaZNmyYsFouorq4Wb775prDZbEIIIT788ENRWFgoTCaTmDdvnrBYLMJms4mlS5eKjz76yOk9qKqqEv/3f/8nPW9sbHS5rwcPHhSpqani0qVLoqmpScyfP7/X9mfRokXi2LFjQgghNm7cKNLT04UQQnz99ddiwYIFoqmpSVgsFuVLg30AAAXpSURBVPHss8+K/fv3CyGE+OSTT6TP0auvviq+++47IYQQVqtVNDc3O637hhvb3Gq1ivT0dFFfXy+mTZsmfdZ27twp1b59+3bp89Vd26W9z2BX9WSf3LltWuNVYwQA0L9/f4waNQoA8Oijj+Lo0aMIDw+XbtocHx+P48eP4/z5861O78jIkSORl5eH/Px8XLp0Cb6+vgCA8PBwDBkyBGq1GpGRkYiJiZG+iVy6dAkA8OOPP+KJJ56QlhUYGNgr+3P8+HE89NBD8PPzg7+/P8aOHXtLmxvPhwwZgsGDByMkJARarRYRERG4fPkyjh49iocffhhBQUFO/fz5558RGxsLX19f9OvXDw8++GCHNd5w/PhxTJgwAQAwZswYBAQEAACOHj2KU6dO4dVXX8XLL7+MH3/8ERcvXsSJEydw7733IigoCD4+Phg3blyryw0PD0dNTQ0+/vhjHDlyBP369XO5r4cOHcKIESOg1+vRr18/xMbG9sr+NDU14cqVK7jnnnsAXP8s3ezee+9Fv379EBQUdMs2v/H5vXr1KkwmEx566CEAgK+vb7uHGvfu3YuXX34Zr732Gmpra2E0GqFSqTB+/HgAwIQJE/Dzzz+3+/7I2S7tfQa7qqf65O5t0xqvOzSkUqm6dXmlpaXYuXMnAGDevHmIi4vDiBEj8P333+P111/HvHnzEB4efssuukqlkp6rVCq0tLR0uT5P9Kej9d7ct9/22+FwQAjRZnvRxmkrGo1Ges1ms7W57taWK4RAfHw8Zs6cecv00tLSVpfR0tKCjIwMANf/oMyYMQOrVq3CkSNHsG/fPhw6dAgGg8HlvvaF/jz33HMubVMAUKvV8PHxkR6318/WHDt2DD/++CP+9re/wc/PDytXrmz1PfhtPd29Xbrz/05P9smd26YtXrdHUFtbi19++QUAUFJSgpiYGNTU1Egj98XFxbjnnnswaNCgVqcDwG233Ybm5mYAwEMPPYRVq1Zh1apViIqKwsWLFxEREYEnn3wSY8eOxa+//upybffffz/27dsnPW9sbOyV/bn77rtRWlqKa9eu4erVq6ioqHC5jwAQExODw4cPo6Ghwamfo0aNQkVFBa5du4bm5mZ8//33AICwsDDpuOm3337b6nLvvvtufPPNNwCA//znP7hy5Yq0vm+//RYWi0Va36VLlxAdHY2ffvoJDQ0NsNvt0nLVarX0HsyYMQP19fVoaWnBww8/jJSUFJw6dcrlvo4fPx6VlZUwmUxobm5GeXl5r+xPQEAA/P39pW+sN9brKn9/f/Tv318KI5vN1uYv0ZqamhAQEAA/Pz+cO3cOJ06cAHA94G7UXFJSIu3p3tCd26W9z2BX9GSf3Llt2uJ1ewR33HEHCgsLsX79egwYMADPP/88oqOjsXbtWmkQNTExEVqtFgaDwWk6ACQkJOCNN95ASEiI0+DqoUOH8M0330Cj0SA4OBhTp05FU1OTS7X94Q9/wEcffYRFixZBrVZj6tSp+N3vftfr+hMYGIjx48fj5ZdfRlhYmNOHuyORkZF4+umnsXLlSqjVagwdOhTz588HAIwYMQIPPvigtOyoqCj4+/sjKSkJ2dnZKC4uxn333dfqcqdNm4Z169YhIyMDd999N/R6PQBg8ODBSElJwd/+9jcIIaDRaDB37lyMHDkS06ZNw9KlSxEcHIxhw4bdsnd2g8lkwvvvvy+99ttv4u155JFHcPDgQSxYsAA+Pj7S7ntv7I/BYJAGJB944AGX+3jDggULsH79euzYsQMajQbp6emIiIhwmm/06NH48ssvsXjxYgwaNAjR0dEAAD8/P1RXVyMjIwP+/v5YuHBhu+uTs13a+wx2RU/3yV3bpk2dGlEgRWpr4Kurrl69KoQQorm5WWRkZIiqqqpuW7Yn9MX+XLx4URqQdJdZs2a5dX3u0BN98sS28bo9Aur9PvjgA5w9exY2mw3x8fG94qJbcnhbf0h5eNE5IiKF87rBYiIi6hwGARGRwjEIiIgUjkFARKRwDAIiIoX7fx/nxahxvNYoAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "results = particles.multiSMC(fk={'boot':fk_gp, 'guided':fk_guided, 'apf': fk_apf}, \n", " N=100, qmc={'smc': False, 'sqmc': True}, nruns=200)\n", "\n", "sb.boxplot(x=['%s-%s'%(r['fk'], r['qmc']) for r in results], y=[r['output'].logLt for r in results]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's discard the bootstrap algorithms to better visualise the results for the other algorithms:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de1BUZ57/8Xd3gygYQeQiKo4G8RKvMWAJwVhsQdWUCW7mDx104o5TmXIVjbU6yYIV1EwYRzLKWElqNqhx3IpjZjZjvMDMrNZSjDpEMILjuMkm+1PQCEJCAGm8gTR9fn8Qe+Pg8UKD3Q2f1z+h+/Cc/vaTI5/znOdcLIZhGIiIiNyF1dMFiIiI91JIiIiIKYWEiIiYUkiIiIgphYSIiJhSSIiIiCk/TxfQ02praz1dgoiITxkxYoTpMo0kRETElEJCRERMKSRERMSUQkJEREwpJB6x5uZmcnNzsdvtni5FROS+3Dq7ac+ePVRUVODn50dkZCQZGRkEBQXhcDjYsWMHlZWVWK1Wli5dyuTJk7u0v3jxIjt37qS1tZXw8HBWr15NYGAgZ8+eZe/evTgcDvz8/FiyZAlTpkxxp1SvUVhYyLlz5ygoKGDJkiWeLkdE5J4s7twF9m9/+xtTpkzBZrPxm9/8BoAXXniBw4cPU1VVRUZGBna7nZ///Ods3rwZq/XOgcu6detYsmQJTzzxBMXFxdTX15Oens6FCxcIDg4mNDSUS5cusWnTJrZv3/5ANXnzKbDNzc1kZmbS3t6Ov78/v/jFLwgODvZ0WSLSz/XaKbDTp0/HZrMBMH78eJqamgCoqalx7fkHBwcTFBREVVVVl/a1tbVMmjQJgGnTpnHy5EkAxo4dS2hoKADR0dG0t7fT3t7uTqleobCwEKfTCYDT6aSgoMDDFYmI3FuPXUxXXFxMYmIiAGPGjKG8vJynn36axsZGqqqqaGhoYNy4cXe0iY6Opry8nPj4eMrKymhsbOyy3pMnTzJ27Fj8/f3v+rlFRUUUFRUBkJubS1hYWE99pR5XVlZGR0cHAB0dHZSVlbFmzRoPVyUiYu6+IZGTk0Nzc3OX99PT04mPjwdg//792Gw25syZA0BycjI1NTVkZWURHh7OhAkTXCOOb1uxYgW7d+9m3759xMXF4ed3ZznV1dXs3buXV1991bS+lJQUUlJSXK8bGhru95U8Zvbs2Rw/fpyOjg5sNhuzZ8/26npFpH+41+Gm+4bE+vXr77n86NGjVFRUsGHDBiwWCwA2m42lS5e6fic7O5uoqKgubUeOHEl2djbQeejp9OnTrmWNjY1s3bqVlStXMnz48PuV6RPS0tIoKSmho6MDq9XK/PnzPV2SiMg9uTUncebMGQ4dOkRmZiYBAQGu99va2mhtbQXg7Nmz2Gw2Ro0a1aX97dNAnU4n+/fvJzU1FYDr16+Tm5vLokWLmDhxojslepWQkBCSkpKwWCwkJSVp0lpEvJ5bZze99NJLOBwOBg8eDEBsbCzLli2jvr6eTZs2YbVaCQ0NZfny5YSHhwOQn59PamoqMTEx/OlPf+LIkSMAzJo1i8WLF2OxWPjwww85ePDgHSOI7OzsB/qj6s1nN0HnGU75+fmsWLFCISEiXuFeh5vcCglv5O0hIeLNtBPTP+kusCLyQL59sacIKCRE5BvNzc2UlJRgGAYlJSW6dYwACgkR+YYu9pS7UUiICAClpaV3XOxZWlrq4YrEGygkRASAhIQE10WvNpuNhIQED1ck3kAh8YjpVuHirdLS0lw34dTFnnKbQuIR09kj4q10safcjULiEdLZI+Lt0tLSiI2N1ShCXBQSj5DOHhFvFxISQlZWlkYR4qKQeIR09oiI+BqFxCOks0dExNcoJB4hnT0iIr6mx55MJ/d3++yRo0eP6uwR6XHvv/8+1dXVbq3jq6++AiAyMrLb64iOjmbx4sVu1SHeQyHxiKWlpXH58mWNIsQrtbW1eboE8TK6VbiIuLzxxhsAZGZmergSeZR0q3AREekWhYSIiJjSnMRD8JaJQdDkoIg8GgqJR0wTgyLiSxQSD6En9tw1MSgivkRzEiIiYkohISIiphQSIiJiSiEhIiKmFBIiImJKISEiIqYUEiIiYkohISIiphQSIiJiyq0rrvfs2UNFRQV+fn5ERkaSkZFBUFAQDoeDHTt2UFlZidVqZenSpUyePLlL+4sXL7Jz505aW1sJDw9n9erVBAYGupY3NDSwZs0aFixYoOcviIh4gFsjiWnTppGXl8fWrVuJioriwIEDABQVFQGQl5dHdnY27733Hk6ns0v77du384Mf/IC8vDxmzZpFQUHBHcv//d//nSeffNKdEkVExA1uhcT06dOx2WwAjB8/nqamJgBqamqYMmUKAMHBwQQFBVFVVdWlfW1tLZMmTQI6A+fkyZOuZR9//DGRkZGMGjXKnRJFRMQNPXaDv+LiYhITEwEYM2YM5eXlPP300zQ2NlJVVUVDQwPjxo27o010dDTl5eXEx8dTVlZGY2MjAK2trRw6dIj169d3GV38vaKiItfIJTc3l7CwsJ76Sr3C398fwOvrlP5J26f8vfuGRE5ODs3NzV3eT09PJz4+HoD9+/djs9mYM2cOAMnJydTU1JCVlUV4eDgTJkxwjTi+bcWKFezevZt9+/YRFxeHn19nOR988AHPPvssAwcOvO8XSElJISUlxfW6oaHhvm08qb29HfD+OqV/0vbZP93r8aX3DYn169ffc/nRo0epqKhgw4YNWCwWAGw2G0uXLnX9TnZ2NlFRUV3ajhw5kuzsbKDz0NPp06cBOH/+PCdPnmTv3r1cv34di8XCgAED+O53v3u/ckVEpAe5dbjpzJkzHDp0iJ/+9KcEBAS43m9ra8MwDAYOHMjZs2ex2Wx3nVuw2+0EBwfjdDrZv38/qampALz++uuu3/nggw8YOHCgAkJExAPcColdu3bhcDjIyckBIDY2lmXLlmG329m0aRNWq5XQ0FBWrVrlapOfn09qaioxMTF89NFHHDlyBIBZs2aRnJzsTjkiItLDLIZhGJ4uoifV1tZ6uoR70pPpxJtp++yf7jUnoSuuRUTElEJCRERMKSRERMSUQkJEREwpJERExJRCQkRETCkkRETElEJCRERMKSRERMSUQkJEREwpJERExJRCQkRETCkkRETElEJCRERMKSRERMSUQkJEREwpJERExJRCQkRETCkkxKc1NzeTm5uL3W73dCkifZJCQnxaYWEh586do6CgwNOliPRJCgnxWc3NzZSUlGAYBiUlJRpNiPQChYT4rMLCQpxOJwBOp1OjCfEqfeVQqEJCfFZpaSkdHR0AdHR0UFpa6uGKRP5PXzkUqpAQn5WQkIDNZgPAZrORkJDg4YpEOvWlQ6EKCfFZaWlpWK2dm7DVamX+/PkerkikU186FKqQEJ8VEhJCUlISFouFpKQkgoODPV2SCNC3DoUqJMSnpaWlERsbq1GEeJWZM2fe8fqpp57yUCXu8/N0ASLuCAkJISsry9NliNyTYRieLqHb3AqJPXv2UFFRgZ+fH5GRkWRkZBAUFITD4WDHjh1UVlZitVpZunQpkydP7tL+4sWL7Ny5k9bWVsLDw1m9ejWBgYEAfPHFF+zYsYObN29isVjYvHkzAwYMcKdcEZFH4vTp0/d87UvcColp06axePFibDYbv/nNbzhw4AAvvPACRUVFAOTl5WG32/n5z3/O5s2bXZOMt23fvp0lS5bwxBNPUFxcTEFBAenp6XR0dPD222+zatUqxowZw9WrV/Hz06BHRHzDzJkzOXHihOu1Lx9ucmtOYvr06a5TEMePH09TUxMANTU1TJkyBYDg4GCCgoKoqqrq0r62tpZJkyYBnYFz8uRJAP72t78xevRoxowZA8Bjjz3WJWBERHxFvz3c9G3FxcUkJiYCMGbMGMrLy3n66adpbGykqqqKhoYGxo0bd0eb6OhoysvLiY+Pp6ysjMbGRgDq6uqwWCxs2rSJlpYWEhMT+cd//Me7fm5RUZFr5JKbm0tYWFhPfaVe4e/vD+D1dUr/pO2zZ/z1r3/t8tpX+/S+IZGTk0Nzc3OX99PT04mPjwdg//792Gw25syZA0BycjI1NTVkZWURHh7OhAkTXCOOb1uxYgW7d+9m3759xMXFuQ4pdXR08Pnnn7N582YCAgJ4/fXXefzxx5k6dWqXdaSkpJCSkuJ63dDQ8IBf3TPa29sB769THq3333+f6upqT5fBpUuXAHjllVc8Wkd0dDSLFy/2aA3umD17NseOHcPpdGK1Wpk9e7ZX/5sfMWKE6bL7hsT69evvufzo0aNUVFSwYcMGLBYL0Hn169KlS12/k52dTVRUVJe2I0eOJDs7G+g89HR7cmfYsGE88cQTDBkyBIAnn3ySCxcu3DUkxHf1xB/Gr776CoDIyEi31uPpP0rV1dVUVX5BSHDXfyePlNE5kmhquOWxEprtdR777J6SlpZGSUkJTqcTm83m06dou3W46cyZMxw6dIif/vSnBAQEuN5va2vDMAwGDhzI2bNnsdlsjBo1qkt7u91OcHAwTqeT/fv3k5qaCnTOdRQUFNDW1oafnx+fffYZzz77rDulSh/V1tbm6RJ6TEhwFP8w58eeLsPjiv/yrqdLANzfibm90xwYGEh+fn631+PpHRi3QmLXrl04HA5ycnIAiI2NZdmyZdjtdjZt2oTVaiU0NJRVq1a52uTn55OamkpMTAwfffQRR44cAWDWrFkkJycDMHjwYJ599lnWrVuHxWLhySef7HJxivi+ntjw33jjDQAyMzPdXpdIT7JarVgsFp+di7jNrZB4++237/p+REQEb7755l2XLV++3PXzvHnzmDdv3l1/75lnnuGZZ55xpzwRkW5zdyemr+zA6LxSEREx1a+uUPOGM0hunz1yey/Dkzx9rFNEvF+/Conq6mq++H//S1TXs3EfGf/OG0Nyq/J/PVcEUNfh0Y8XER/Rr0ICIMoGPw7WUbZ37U5PlyAiPkB/LUVExJRCQkRETCkkRETElEJCRERMKSRERMSUQkJEREwpJERExJRCQkRETCkkRETEVL+74lp6hjfcBwu8515Yug+Wd/GG7dNbtk1wb/tUSEi3VFdX8/n5CzB4mGcL6eh8sMvnX7Z4roZrjZ77bLmr6upqzp07x6BBgzxWw+1HFdfU1HisBoCbN2+61V4hId03eBjGjOc8XYXHWc78wdMlyF0MGjSI2NhYT5fhcefOnXOrveYkRETElEJCRERM6XCTiPQ59fX13Lx50+1DLX3BzZs3qa+v73Z7hYR0S319PVy7ruPxANcaqa9v9XQVIr2iX4VEfX09Nx164A5AnQMGubF3IeLNIiIiuHXrliau6Zy4joiI6Hb7fhUS0nMiIiJocrbo7CY6z26KiBji1jrq6+u5dvUmxX95t4eq8l3N9jocTs+duip36lchERERwa2rV/T4UjpHUwPc2LsQkf6hX4WEiLeKiIjAz3qLf5jzY0+X4nHFf3mX0LABni5DvqFdahERMaWRhHTftUbPn910097530HBnqvhWiPg3pyE9DxPnwLb1tYGQEBAgMdqAN2WQzwkOjra0yUAcOlSMwCjh3vyj/QQr+kP6eQN/z9u3+Bv1KhRHq7Evf5QSEi3eMsdT2/fYTMzM9PDlYg38Ybts69sm5qTEBERU26NJPbs2UNFRQV+fn5ERkaSkZFBUFAQDoeDHTt2UFlZidVqZenSpUyePLlL+4sXL7Jz505aW1sJDw9n9erVBAYG4nA4yM/P58KFCzidTp555hm+973vuVOqiIh0g1sjiWnTppGXl8fWrVuJioriwIEDABQVFQGQl5dHdnY27733Hk5n16uct2/fzg9+8APy8vKYNWsWBQUFAJSVleFwOMjLyyM3N5eioiK37j0iIiLd41ZITJ8+HZvNBsD48eNpamoCOh+yMWXKFACCg4MJCgqiqqqqS/va2lomTZoEdAbOyZMnXctaW1vp6Ojg1q1b+Pn5ERgY6E6pIiLSDT02cV1cXExiYiIAY8aMoby8nKeffprGxkaqqqpoaGhg3Lhxd7SJjo6mvLyc+Ph4ysrKaGzsfMLX7NmzKS8vZ9myZdy6dYsf/vCHDB48+K6fW1RU5Bq55ObmEhYWZlqjv78/t3riy/YR/v7+9+wvX+Dv7w/QR76Hts7btG16j/uGRE5ODs3NzV3eT09PJz4+HoD9+/djs9mYM2cOAMnJydTU1JCVlUV4eDgTJkxwjTi+bcWKFezevZt9+/YRFxeHn19nOefPn8dqtbJ9+3auX7/Ohg0bmDp1KpGRkV3WkZKSQkpKiut1Q0OD6Xe5/ThB6dTe3n7P/vIFt/+f9pXvIZ20bT5aI0aMMF1235BYv379PZcfPXqUiooKNmzYgMXS+bxhm83G0qVLXb+TnZ1NVFRUl7YjR44kOzsb6Dz0dPr0aQBKSkqYMWMGfn5+BAcHM2HCBCorK+8aEiIi0nvcmpM4c+YMhw4dIjMz846rCtva2mht7by//tmzZ7HZbHe9oMRu77xa1ul0sn//flJTU4HO4dknn3yCYRi0trZy7tw5Ro4c6U6pIiLSDW7NSezatQuHw0FOTg4AsbGxLFu2DLvdzqZNm7BarYSGhrJq1SpXm/z8fFJTU4mJieGjjz7iyJEjAMyaNYvk5GQAvvvd7/Jv//Zv/OQnP8EwDJKTk/nOd77jTqkiItINboXE22+/fdf3IyIiePPNN++6bPny5a6f582bx7x587r8zsCBA1m7dq07pYmISA/QFdciImJKISEiIqYUEiIiYkp3gRXxEs32Oo8/4/ra9c4LWgcHDfNYDc32OkLDdKKKt1BIiE9zOBzU1tZit9sJDvbgg4fc5A3PPwC4dqPzAjBPPj40NOw7XtMfopAQD3r//feprq52ax0XL16ko6OD1157jeHDh3d7PdHR0R59BoE3PP8A+s4zEHqCu9vnF198QVtbG5s2bXLdTaI7PL1tak5CfJbD4aCjowOAlpYWHA6HhysS+T9OpxPDMHzithz30u9GEnUd8K69623LH5XGzr9pDOt6K6tHqq4DPH3U1929oz179rhGElarldGjR7NkyZIeqk76O3e2z+bmZtdo7MaNGyxfvtxnD4f2q5DwhuOc7d8893bA6NEereM7eEd/uKO0tNQ1kujo6KC0tFQhIV6hsLDQ9Qwdp9NJQUGBz26b/SokvOG4r4759pyEhASOHz9OR0cHNpuNhIQET5ckAvStHRjNSYjPSktLw2rt3IStVivz58/3cEUinRISElyPR/D1HRiFhPiskJAQkpKSsFgsJCUl+ewxX+l7+tIOjEJCfFpaWhqxsbE+/Y9Q+p6+tAPTr+YkpO8JCQkhKyvL02WIdJGWlsbly5d9fgdGISEi0gv6yg6MDjeJiIgphYSIiJhSSIiIiCmFhIiImFJIiIiIKYWEiIiYUkiIiIgphYSIiJhSSIiIiCmFhIiImFJIiIiIKYWEiIiYUkiIiIgpt+4C+7vf/Y7y8nIsFgvBwcFkZGQQGhoKwIEDByguLsZqtfKjH/2IGTNmdGl/7do1tm3bxtdff014eDhr1qxh8ODBD9xeRER6l1sjifnz57N161a2bNnCzJkz2bdvHwA1NTWcOHGCX/7yl7z66qvs2rXL9VDwbzt48CBTp07lrbfeYurUqRw8ePCh2ouISO9yKyQCAwNdP7e1tWGxWAA4deoUiYmJ+Pv7ExERwfDhwzl//nyX9qdOnWLu3LkAzJ07l1OnTj1UexER6V1uP3Tot7/9LcePHycwMJCNGzcC0NTURGxsrOt3QkNDaWpq6tLWbrczdOhQAIYOHUpLS8tDtQcoKiqiqKgIgNzcXMLCwtz9Sr3K398fwOvrlP5J26f8vfuGRE5ODs3NzV3eT09PJz4+nkWLFrFo0SIOHDjA4cOHWbhwIYZhuFXUw7RPSUkhJSXF9bqhocGtz+5t7e3tgPfXKf2Tts/+acSIEabL7hsS69evf6APSUpKIjc3l4ULFzJs2DAaGxtdy5qamlwT2t8WHBzMlStXGDp0KFeuXGHIkCEAD9xeRER6l1tzEnV1da6fy8vLXWkUFxfHiRMnaG9vp76+nrq6OsaNG9elfVxcHMeOHQPg2LFjxMfHP1R7ERHpXW7NSezdu5e6ujosFgthYWEsW7YMgOjoaBISEli7di1Wq5UXX3wRq7Uzj/Lz80lNTSUmJobnn3+ebdu2UVxcTFhYGGvXrr1vexEReXQshrsTCF6mtrbW0yXc0xtvvAFAZmamhysR6UrbZ/90rzkJ7Z6LiIgphYSIiJhSSIiIiCmFhIiImFJIiIiIKYWEiIiYUkiIiIgphYSIiJhSSIiIiCmFhIiImFJIiIiIKYWEiIiYUkiIiIgphYSIiJhSSIiIiCmFhIiImFJIiIiIKYWEiIiYUkiIiIgphYSIiJhSSIiIiCmFhIiImFJIiIiIKYWEiIiYUkiIiIgphYSIiJjy83QBItIz3n//faqrq91ax6VLlwB44403ur2O6OhoFi9e7FYd4j0UEiLiEhAQ4OkSxMsoJET6CO29S29wKyR+97vfUV5ejsViITg4mIyMDEJDQwE4cOAAxcXFWK1WfvSjHzFjxowu7a9du8a2bdv4+uuvCQ8PZ82aNQwePJizZ8+yd+9eHA4Hfn5+LFmyhClTprhTqoiIdINbE9fz589n69atbNmyhZkzZ7Jv3z4AampqOHHiBL/85S959dVX2bVrF06ns0v7gwcPMnXqVN566y2mTp3KwYMHAXjsscfIzMwkLy+PlStX8vbbb7tTpoiIdJNbIREYGOj6ua2tDYvFAsCpU6dITEzE39+fiIgIhg8fzvnz57u0P3XqFHPnzgVg7ty5nDp1CoCxY8e6RiTR0dG0t7fT3t7uTqkiItINbs9J/Pa3v+X48eMEBgayceNGAJqamoiNjXX9TmhoKE1NTV3a2u12hg4dCsDQoUNpaWnp8jsnT55k7Nix+Pv73/Xzi4qKKCoqAiA3N5ewsDB3v1Kvuv09vL1OERF4gJDIycmhubm5y/vp6enEx8ezaNEiFi1axIEDBzh8+DALFy7EMIweKa66upq9e/fy6quvmv5OSkoKKSkprtcNDQ098tm95faIyNvrFJH+Y8SIEabL7hsS69evf6APSUpKIjc3l4ULFzJs2DAaGxtdy5qamlyHj74tODiYK1euMHToUK5cucKQIUNcyxobG9m6dSsrV65k+PDhD1SDiIj0LLcON9XV1REVFQVAeXm5K43i4uJ46623eO6557hy5Qp1dXWMGzeuS/u4uDiOHTvG888/z7Fjx4iPjwfg+vXr5ObmsmjRIiZOnOhOiT3KWy5WAl2wJCKPhlshsXfvXurq6rBYLISFhbFs2TKg8w9YQkICa9euxWq18uKLL2K1ds6R5+fnk5qaSkxMDM8//zzbtm2juLiYsLAw1q5dC8Dhw4f58ssv+fDDD/nwww8ByM7OJjg42J1yvYIuVhIRX2IxemoCwUvU1tZ6ugQREZ9yrzkJ3eBPRFyam5vJzc3Fbrd7uhTxEgoJEXEpLCzk3LlzFBQUeLoU8RIKCREBOkcRJSUlGIZBSUmJRhMCKCRE5BuFhYWu2+c4nU6NJgRQSIjIN0pLS+no6ACgo6OD0tJSD1ck3kAhISIAJCQkYLPZALDZbCQkJHi4IvEGCgkRASAtLc11PZPVamX+/Pkerki8gUJCRAAICQkhKSkJi8VCUlJSn7h4VdynJ9OJiEtaWhqXL1/WKEJcdMW1iEg/pyuuRUSkWxQSIiJiSiEhIiKmFBIiImKqz01ci4hIz9FIwgOysrI8XUKfov7sWerPntMX+lIhISIiphQSIiJiSiHhASkpKZ4uoU9Rf/Ys9WfP6Qt9qYlrERExpZGEiIiYUkiIiIgphcRDqqys5Ne//vVdl61cuZKWlpYHXtfRo0fZtWtXT5XmM9SHnnX58mVeeeUV/vVf/5Uvv/zS0+X4pP7Uh7pV+EOKiYkhJibG02X4NPWhZ506dYr4+HgWLlzo6VJ8Vn/qQ4UEsG/fPkpKShg2bBiPPfYYjz/+OKdPn2bJkiXExMTQ0tLCunXr+NWvfsWnn35KYWEhWVlZXL16lTfffJOWlhbGjRvHt88BOH78OP/5n/+Jw+EgNjaWH//4x1itVv785z9z8OBBQkJCiIqKwt/fv0s9TqeTd955h6qqKgCSk5N57rnneO211xgzZgwXLlygpaWFlStXcvDgQS5dukRiYiLp6ekAHDt2jMLCQiwWC6NHj+all15SH37Th1VVVbzzzjsMGDCAiRMncubMGfLy8jh69Cgff/wxTqeT6upq0tLScDgcHD9+HH9/f9atW8fgwYP58ssv2blzJy0tLVitVtasWcPw4cN7vX/v5xe/+AWNjY20t7czb948UlJSWLJkCampqXz66acEBQXxL//yL5w/f54//vGPWK1WPvvsMzZu3Ohahy9udz2pN/vQp7c7o587f/688fLLLxttbW3GjRs3jJdeesk4dOiQsXHjRuP8+fOGYRiG3W43MjIyDMMwjE8++cTYvHmzYRiGsWvXLuP3v/+9YRiGUVFRYSxYsMCw2+1GdXW1sXnzZqO9vd0wDMPYuXOncfToUaOpqclYvny5Ybfbjfb2diM7O9t49913u9RUWVlpvP76667X165dMwzDMDZu3Gjs2bPHMAzD+OMf/2gsW7bMaGpqMm7dumX88z//s9HS0mJcunTJWL16tWG32w3DMIyrV6/2RrfdwZf68Cc/+Ynx6aefGoZhGO+9956xdu1awzAM489//rOxatUq48aNG4bdbjf+6Z/+yThy5IhhGIaxe/du4w9/+INhGIaxbt064+TJk4ZhGEZbW5vR2tragz3Zfbf/P7e1tRlr1641WlpajAULFhjHjx83DMMwfv/737v66T/+4z+MQ4cOdVmHr213Pa03+9CXt7t+P5L4/PPPiY+PZ8CAAQA89dRTD9z2s88+4+WXXwZg5syZBAUFAfDJJ59w4cIF1q1bB8CtW7cYMmQI586dY/LkyQwZMgTofPB8XV1dl/VGRERQX1/Pr3/9a2bOnMm0adNcy+Li4gAYPXo0o0aNYujQoQBERkwF2gcAAANZSURBVEbS2NjIZ599xuzZs12fMXjw4Ifqj+7wlT68ceMG169f54knngDgmWee4cyZM642kydPZtCgQQwaNIjAwMA7+vrSpUvcvHmTpqYmZs2aBeD6vt7gT3/6E6dOnQKgoaGBuro6LBYLiYmJAMyZM4etW7fecx2+tt31tN7qQ1/f7vp9SBgml4nYbDbXsvb2dtP2FovlruucO3cuixcvvuP9jz/++K7rcDqdZGZmAp3/GL///e+zZcsWzpw5w+HDhzlx4gQZGRkArkMrFovljsMsFouFjo4ODMO4a029yVf68Ic//OE9++bb/Wm1WvHz83P9fLtvvdGnn37Kf//3f/Ozn/2MgIAAXnvttbv2999/d1/f7npSb/ahr293/f7spokTJ1JRUcGtW7dobW3l9OnTAISHh7uOK5aVld217aRJk/jLX/4CwF//+leuX78OwNSpUykrK8NutwNw7do1vv76a2JjY/mf//kfrl69isPhcK3XarWyZcsWtmzZwve//31aWlpwOp3Mnj2b9PR0Lly48MDfZ+rUqZSWlnL16lXXZ/c2X+nDoKAgAgMD+fzzzwFcn/ugAgMDGTZsmCuo2tvbaWtre6h19IYbN24QFBREQEAAly9f5ty5c0Bn0N7un5KSEiZOnHhHO1/f7npSb/ahr293/X4kMW7cOJ566ileeeUVwsPDiYmJITAwkLS0NLZt28bx48eZMmXKXdsuWLCAN998k8zMTCZNmkRYWBgAo0aNIj09nZ/97GcYhoHNZuPFF19k/PjxLFiwgOzsbEJCQhg7dixOp7PLepuamnjnnXdcy/5+b/peoqOj+d73vsdrr72G1WplzJgxrFy5shs98+B8qQ8zMjJcE4jTp09/6O+6atUqduzYwQcffIDNZmPt2rVERkY+9Hp60owZM/iv//ovXn75ZUaMGEFsbCwAAQEBVFdXk5mZSWBgIGvWrLnnenxtu+tJvd2HPr3dPbLZDy928+ZNwzAMo7W11cjMzDQqKys9XJHv8cU+/Oqrr1wTiH3RCy+84OkSfF5v9KGvbXf9fiQBsH37dmpqamhvb2fu3Lk8/vjjni7J56gPRfom3eBPRERM9fuJaxERMaeQEBERUwoJERExpZAQERFTCgkRETH1/wEAaiq3yBVP9gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "res_noboot = [r for r in results if r['fk']!='boot']\n", "sb.boxplot(x=['%s-%s'%(r['fk'], r['qmc']) for r in res_noboot], y=[r['output'].logLt for r in res_noboot]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "VoilĂ !" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.3" } }, "nbformat": 4, "nbformat_minor": 1 }