Skip to content

Instantly share code, notes, and snippets.

@simeoncarstens

simeoncarstens/smc.ipynb

Last active January 27, 2020 14:21
Show Gist options
  • Save simeoncarstens/3d3cc9bb340edeecf21bb58f419da860 to your computer and use it in GitHub Desktop.
Save simeoncarstens/3d3cc9bb340edeecf21bb58f419da860 to your computer and use it in GitHub Desktop.
Download ZIP
Importance Sampling, Sequential Monte Carlo, Particle Filtering and stuff
Display the source blob
Display the rendered blob
Raw
smc.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Importance Sampling, Sequential Importance Sampling and Particle Filters\n",
"Matthias, Yves, Leo and I, we recently discussed above sampling algorithms in the context of `monad-bayes`. To get a better understanding, I decided to implement them in Python."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Importance sampling\n",
"... is a very general sampling method. For once, it has nothing to do with Markov chains and not even with Monte Carlo. \n",
"Let's say you want to calculate the expectation of a function $f(x)$ w.r.t. a probability distribution $p(x)$. It is given by the following integral:\n",
"\\begin{equation}\n",
" \\langle f \\rangle_{p} = \\int \\mathrm d x \\ f(x) p(x)\n",
"\\end{equation}\n",
"Now in general, you might not be able to calculate this analytically, so you resort to numerics. If you could sample correctly from $p(x)$, you could draw $N$ samples $x_1, \\ldots, x_N$ and approximate the integral by a sum:\n",
"\\begin{equation}\n",
" \\langle f \\rangle_{p} \\approx \\frac{1}{N} \\sum_{i=1}^N f(x_i) p(x_i)\n",
"\\end{equation}\n",
"The world unfortunately sometimes is kind of a bad place and sampling from $p(x)$ might not be possible. But you might be able to sample from a different distribution $q(x)$. Then, doing what physicists have been doing all their lifes, you insert a 1 and write:\n",
"\\begin{align}\n",
" \\langle f \\rangle_{p} &= \\int \\mathrm d x \\ f(x) \\frac{p(x)}{q(x)}q(x) \\\\\n",
" &= \\int \\mathrm d x \\ f(x) w(x) q(x) \\\\\n",
" &= \\langle f\\cdot w \\rangle_q\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, we introduced the *importance weight* $w(x)=\\frac{p(x)}{q(x)}$ and we see that we're now calculating the expectation value of a different functon w.r.t. a different distribution. That expectation value can now be again approximated, supposed we draw $N$ samples $x_1, \\ldots, x_N$ from $q$ (confusing, right?):\n",
"\\begin{equation}\n",
" \\langle f \\rangle_p = \\langle f\\cdot w \\rangle_q \\approx \\frac{1}{N}\\sum_{i=1}^N f(x_i) w(x_i)\n",
"\\end{equation}\n",
"TODO: I'm not clear why we have to normalize the importants weights (which I haven't done here, but later in the code I do)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Voili, voilà—this is called *importance sampling*, because instead from $p$, we sample from an *importance distribution* $q$. Now you see a problem here: if the distributions $p$ and $q$ don't match up well (and in general, they won't), you will draw many samples from $q$ which have low importance weights $w(x_i)$ and thus contribute only little to the approximation, meaning you will need a lot of samples and still might get a bad estimate."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Implementing importance sampling\n",
"Let's first import useful stuff and define a few PDFs:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.special import logsumexp\n",
"choice = np.random.choice\n",
"\n",
"\n",
"class Pdf(object):\n",
"\n",
" def __call__(self, x):\n",
" \"\"\"log-probability\"\"\"\n",
" pass\n",
"\n",
" def sample(self, n):\n",
" pass\n",
"\n",
"\n",
"\n",
"class Gaussian(Pdf):\n",
"\n",
" def __init__(self, mu=0, sigma=1):\n",
"\n",
" self.mu = mu\n",
" self.sigma = sigma\n",
"\n",
" def __call__(self, x):\n",
"\n",
" return -0.5 * (x - self.mu) ** 2 / self.sigma / self.sigma\n",
"\n",
" def sample(self, n):\n",
"\n",
" return np.random.normal(self.mu, self.sigma, n)\n",
"\n",
"\n",
"class Uniform(Pdf):\n",
"\n",
" def __init__(self, low, high):\n",
"\n",
" self.low = low\n",
" self.high = high\n",
"\n",
" def __call__(self, x):\n",
"\n",
" return np.repeat(-np.log(self.high - self.low), len(x))\n",
"\n",
" def sample(self, n):\n",
"\n",
" return np.random.uniform(self.low, self.high, n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we whip up a class implementing Importance Sampling. Its `sample` method will return the samples from the importance distribution and the corresponding importance weights. Also, as you might have noticed, we work in log-space for numeric stability."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"class ImportanceSampler(object):\n",
"\n",
" def __init__(self, posterior):\n",
"\n",
" self.posterior = posterior\n",
"\n",
" def sample(self, imp_dist, n_samples):\n",
"\n",
" imp_samples = imp_dist.sample(n_samples)\n",
" imp_weights = self.calculate_weights(imp_dist, imp_samples)\n",
" # doing stuff in log-space to avoid underflows\n",
" normalized_imp_weights = imp_weights - logsumexp(imp_weights)\n",
"\n",
" return imp_samples, normalized_imp_weights\n",
"\n",
" def calculate_weights(self, imp_dist, imp_samples):\n",
"\n",
" return self.posterior(imp_samples) - imp_dist(imp_samples)\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's try this out. We want to sample from a Gaussian distribution, but we pretend we cannot, so we sample from a uniform importance distribution instead:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"n_samples = 200000\n",
"target = Gaussian()\n",
"imp = Uniform(-10, 10)\n",
"sampler = ImportanceSampler(target)\n",
"biased_samples, logws = sampler.sample(imp, n_samples)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now I've been talking a lot about approximating expectation values, but for personal reasons, I'm actually more interested in the actual samples. So how do you get samples from the target distribution from the importance samples and their weights? Well, you just draw from the set of importance samples, but with probabilities given by the importance weights:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"samples = choice(biased_samples, n_samples, p=np.exp(logws))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Gesagt, getan, so now we can plot a histogram of the actual (unbiased) samples:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAPu0lEQVR4nO3df6zddX3H8edrrTgzNVTbdaRtVqL9p7qt6k1pwpIxWUoBs2JiDCyTzhFrYkkgYZlV/6hRSWoWdSNTlqqNJWNWIhoaqasdITH+AfaCFSjIuMEy2hR6tSguJpq69/44nzvO6u29p/fX97b3+UhOzve8v5/v97zPN9DX/f46J1WFJGlh+52uG5Akdc8wkCQZBpIkw0CShGEgSQIWd93AVC1durRWr17ddRuSdF555JFHflJVy86sn7dhsHr1aoaHh7tuQ5LOK0meG6/uYSJJkmEgSRogDJKsSvJgkieTHElyS6t/PMnxJIfb45q+ZT6SZCTJ00mu6qtvarWRJNv76pcmebjVv5bkopn+oJKksxtkz+A0cFtVrQU2ANuSrG3zPldV69pjP0Cbdz3wFmAT8IUki5IsAj4PXA2sBW7oW8+n27reDLwE3DRDn0+SNIBJw6CqTlTVo236F8BTwIoJFtkM7K2qX1XVj4ERYH17jFTVs1X1a2AvsDlJgHcCX2/L7wGum+oHkiSdu3M6Z5BkNfA24OFWujnJY0l2J1nSaiuA5/sWO9ZqZ6u/EfhZVZ0+oz7e+29NMpxkeHR09FxalyRNYOAwSPJa4F7g1qp6GbgTeBOwDjgBfGZWOuxTVbuqaqiqhpYt+63LZCVJUzTQfQZJXkUvCO6uqm8AVNWLffO/CHyrvTwOrOpbfGWrcZb6T4GLkyxuewf94yVJc2CQq4kCfBl4qqo+21e/pG/Yu4En2vQ+4Pokr05yKbAG+D5wCFjTrhy6iN5J5n3V+0GFB4H3tOW3APdN72NJks7FIHsGlwPvAx5PcrjVPkrvaqB1QAFHgQ8CVNWRJPcAT9K7EmlbVf0GIMnNwAFgEbC7qo609X0Y2JvkU8AP6IWPdF5avf3+gcYd3XntLHciDW7SMKiq7wEZZ9b+CZa5Hbh9nPr+8ZarqmfpXW0kSeqAdyBLkgwDSZJhIEnCMJAkcR7/noE01wa9Skg6H7lnIEkyDCRJhoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJ+OM2Umc/WjPo+x7dee0sdyK5ZyBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSGCAMkqxK8mCSJ5McSXJLq78hycEkz7TnJa2eJHckGUnyWJK3961rSxv/TJItffV3JHm8LXNHkszGh5UkjW+QPYPTwG1VtRbYAGxLshbYDjxQVWuAB9prgKuBNe2xFbgTeuEB7AAuA9YDO8YCpI35QN9ym6b/0SRJg5o0DKrqRFU92qZ/ATwFrAA2A3vasD3AdW16M3BX9TwEXJzkEuAq4GBVnaqql4CDwKY27/VV9VBVFXBX37okSXPgnM4ZJFkNvA14GFheVSfarBeA5W16BfB832LHWm2i+rFx6uO9/9Ykw0mGR0dHz6V1SdIEBg6DJK8F7gVuraqX++e1v+hrhnv7LVW1q6qGqmpo2bJls/12krRgDBQGSV5FLwjurqpvtPKL7RAP7flkqx8HVvUtvrLVJqqvHKcuSZojg1xNFODLwFNV9dm+WfuAsSuCtgD39dVvbFcVbQB+3g4nHQA2JlnSThxvBA60eS8n2dDe68a+dUmS5sAgP3t5OfA+4PEkh1vto8BO4J4kNwHPAe9t8/YD1wAjwC+B9wNU1akknwQOtXGfqKpTbfpDwFeA1wDfbg9J0hyZNAyq6nvA2a77v3Kc8QVsO8u6dgO7x6kPA2+drBdJ0uzwDmRJkmEgSTIMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJDHY11FI6tDq7fcPPPbozmtnsRNdyAwDXbDO5R9RaaHzMJEkyTCQJBkGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CSxABhkGR3kpNJnuirfTzJ8SSH2+OavnkfSTKS5OkkV/XVN7XaSJLtffVLkzzc6l9LctFMfkBJ0uQG2TP4CrBpnPrnqmpde+wHSLIWuB54S1vmC0kWJVkEfB64GlgL3NDGAny6revNwEvATdP5QJKkczdpGFTVd4FTA65vM7C3qn5VVT8GRoD17TFSVc9W1a+BvcDmJAHeCXy9Lb8HuO4cP4MkaZqmc87g5iSPtcNIS1ptBfB835hjrXa2+huBn1XV6TPq40qyNclwkuHR0dFptC5J6jfVMLgTeBOwDjgBfGbGOppAVe2qqqGqGlq2bNlcvKUkLQiLp7JQVb04Np3ki8C32svjwKq+oStbjbPUfwpcnGRx2zvoHy9JmiNT2jNIcknfy3cDY1ca7QOuT/LqJJcCa4DvA4eANe3KoYvonWTeV1UFPAi8py2/BbhvKj1JkqZu0j2DJF8FrgCWJjkG7ACuSLIOKOAo8EGAqjqS5B7gSeA0sK2qftPWczNwAFgE7K6qI+0tPgzsTfIp4AfAl2fs00mSBjJpGFTVDeOUz/oPdlXdDtw+Tn0/sH+c+rP0rjaSJHXEO5AlSYaBJMkwkCRhGEiSMAwkSUzxpjOpS6u33991C9IFxz0DSZJhIEkyDCRJeM5AuqAMej7l6M5rZ7kTnW/cM5AkGQaSJMNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSQwQBkl2JzmZ5Im+2huSHEzyTHte0upJckeSkSSPJXl73zJb2vhnkmzpq78jyeNtmTuSZKY/pCRpYoPsGXwF2HRGbTvwQFWtAR5orwGuBta0x1bgTuiFB7ADuAxYD+wYC5A25gN9y535XpKkWTZpGFTVd4FTZ5Q3A3va9B7gur76XdXzEHBxkkuAq4CDVXWqql4CDgKb2rzXV9VDVVXAXX3rkiTNkameM1heVSfa9AvA8ja9Ani+b9yxVpuofmycuiRpDk37BHL7i75moJdJJdmaZDjJ8Ojo6Fy8pSQtCFMNgxfbIR7a88lWPw6s6hu3stUmqq8cpz6uqtpVVUNVNbRs2bIpti5JOtNUw2AfMHZF0Bbgvr76je2qog3Az9vhpAPAxiRL2onjjcCBNu/lJBvaVUQ39q1LkjRHFk82IMlXgSuApUmO0bsqaCdwT5KbgOeA97bh+4FrgBHgl8D7AarqVJJPAofauE9U1dhJ6Q/Ru2LpNcC320MLzOrt93fdgrSgpXfI//wzNDRUw8PDXbehGWIYzE9Hd17bdQuaYUkeqaqhM+vegSxJMgwkSYaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAlY3HUDurCt3n5/1y1IGoB7BpIkw0CS5GEiSRMY9DDf0Z3XznInmm3uGUiSDANJkmEgSWKaYZDkaJLHkxxOMtxqb0hyMMkz7XlJqyfJHUlGkjyW5O1969nSxj+TZMv0PpIk6VzNxJ7Bn1fVuqoaaq+3Aw9U1RrggfYa4GpgTXtsBe6EXngAO4DLgPXAjrEAkSTNjdk4TLQZ2NOm9wDX9dXvqp6HgIuTXAJcBRysqlNV9RJwENg0C31Jks5iumFQwHeSPJJka6str6oTbfoFYHmbXgE837fssVY7W/23JNmaZDjJ8Ojo6DRblySNme59Bn9aVceT/D5wMMmP+mdWVSWpab5H//p2AbsAhoaGZmy9krTQTWvPoKqOt+eTwDfpHfN/sR3+oT2fbMOPA6v6Fl/ZamerS5LmyJTDIMnvJXnd2DSwEXgC2AeMXRG0BbivTe8DbmxXFW0Aft4OJx0ANiZZ0k4cb2w1SdIcmc5houXAN5OMreffqurfkxwC7klyE/Ac8N42fj9wDTAC/BJ4P0BVnUrySeBQG/eJqjo1jb4kSedoymFQVc8CfzJO/afAlePUC9h2lnXtBnZPtRdJ0vR4B7IkyTCQJBkGkiQMA0kShoEkCX/pTFPkD91LFxb3DCRJ7hlImj5/K/n8556BJMkwkCQZBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJv7VUZ/B3CjSb/HbT+cs9A0mSYSBJMgwkSRgGkiQMA0kShoEkCS8tXTC8ZFTSRAwDSfOO9yPMPQ8TSZIMA0mSh4nOe54LkDQT5k0YJNkE/BOwCPhSVe3suCVJ89y5/DHk+YWJzYvDREkWAZ8HrgbWAjckWdttV5K0cMyXPYP1wEhVPQuQZC+wGXiy065mmId0pO7M9P9/F9qexnwJgxXA832vjwGXnTkoyVZga3v530menoPeJrMU+EnXTcwDbodXuC16LujtkE+f0/D5tC3+cLzifAmDgVTVLmBX1330SzJcVUNd99E1t8Mr3BY9bodXnA/bYl6cMwCOA6v6Xq9sNUnSHJgvYXAIWJPk0iQXAdcD+zruSZIWjHlxmKiqTie5GThA79LS3VV1pOO2BjWvDlt1yO3wCrdFj9vhFfN+W6Squu5BktSx+XKYSJLUIcNAkmQYzJQktyWpJEu77qUrSf4hyY+SPJbkm0ku7rqnuZRkU5Knk4wk2d51P11JsirJg0meTHIkyS1d99SlJIuS/CDJt7ruZSKGwQxIsgrYCPxX17107CDw1qr6Y+A/gY903M+c8StV/p/TwG1VtRbYAGxbwNsC4Bbgqa6bmIxhMDM+B/w9sKDPxlfVd6rqdHv5EL37RRaK//tKlar6NTD2lSoLTlWdqKpH2/Qv6P1DuKLbrrqRZCVwLfClrnuZjGEwTUk2A8er6odd9zLP/C3w7a6bmEPjfaXKgvwHsF+S1cDbgIe77aQz/0jvD8X/6bqRycyL+wzmuyT/AfzBOLM+BnyU3iGiBWGibVFV97UxH6N3qODuuexN80uS1wL3ArdW1ctd9zPXkrwLOFlVjyS5out+JmMYDKCq/mK8epI/Ai4FfpgEeodFHk2yvqpemMMW58zZtsWYJH8DvAu4shbWTSx+pUqfJK+iFwR3V9U3uu6nI5cDf5nkGuB3gdcn+deq+uuO+xqXN53NoCRHgaGqmi/fTjin2g8UfRb4s6oa7bqfuZRkMb2T5lfSC4FDwF+dR3fSz5j0/jLaA5yqqlu77mc+aHsGf1dV7+q6l7PxnIFm0j8DrwMOJjmc5F+6bmiutBPnY1+p8hRwz0IMguZy4H3AO9t/B4fbX8eax9wzkCS5ZyBJMgwkSRgGkiQMA0kShoEkCcNAkoRhIEkC/hdVVliRb+I2xgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"ax.hist(samples, bins=30)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So that's good: looks pretty normal to me."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sequential Importance Sampling\n",
"Now imagine we are dealing with a problem in which data is sequential, like a time series or a linear sequence of nucleotides in DNA. Let's assume that after $n$ (time) steps, you have $n$ observations $y_1, \\ldots, y_n \\equiv y_{1:n}$ and you want to know which parameters $x_1, \\ldots, x_n\\equiv x_{1:n}$ at steps $1, \\ldots, n$ gave rise to these observations. Then you would have to draw samples from the posterior distribution $p(x_{1:n}|y_{1:n})$ of \"paths\" $x_{1:n}$. Say you somehow managed to do that and proceed to step $n+1$. Then you will have to sample from the posterior distributions of paths $x_{1:n+1}$. This gets pretty unwieldy for long sequences. \n",
"But if you can assume that the sequence of $x$'ses is a Markov chain, namely, $p(x_{n+1}|x_{1:n}) = p(x_{n+1}|x_n)$, then you can recycle samples from previous time steps.\n",
"When using Importance Sampling to sample from these distributions, the importance weights can be calculated recursively. A nice special case is obtained if the importance distribution $q$ is just the prior distribution $p(x_{1:n})$. In that case, the importance weight of the $i$-th sample is given by\n",
"\\begin{equation}\n",
" w(x^i_n) \\propto w(x^i_{n-1}) p(y_n|x^i_n)\n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A problem with SIS is now that, after already a few steps, only very few weights will have values significantly different from zero. TODO: why exactly is that? I suspect that somehow the weights keep the memory of the initial stage when there was only very little data and the probabilities $p(y|x)$ are very low becaues of bad estimates of the $x$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So let's implement SIS. First, we write a class representing a \"sequential prior distribution\", which encodes the Markov property. Basically, it draws new samples based only on the most recent batch of samples:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"class SequentialPrior(Pdf):\n",
"\n",
" def __init__(self, distribution, first_state):\n",
"\n",
" self.distribution = distribution\n",
" self.samples = first_state\n",
"\n",
" def sample(self, _):\n",
"\n",
" self.samples = self.distribution.sample(self.samples)\n",
"\n",
" return self.samples"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To instantiate a `SequentialPrior` object, we require some conditional distribution (`distribution`) and an initial batch of samples (`first_state`). The `sample` method then draws samples from the distribution conditioned on the previous batch of samples."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next up is the actual Sequential Importance Sampler. It derives from `ImportanceSampler`, but saves the current weights and modifies the `calculate_weights` method to calculate weights recursively:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"class SequentialImportanceSampler(ImportanceSampler):\n",
"\n",
" def __init__(self, posterior):\n",
"\n",
" super(SequentialImportanceSampler, self).__init__(posterior)\n",
" self.weights = None\n",
"\n",
" def sample(self, n_samples):\n",
"\n",
" (imp_samples, weights) = super(\n",
" SequentialImportanceSampler, self).sample(self.posterior.prior,\n",
" n_samples)\n",
" # store current weights to make them available to recursion\n",
" self.weights = weights\n",
"\n",
" return imp_samples, weights\n",
"\n",
" def calculate_weights(self, _, imp_samples):\n",
"\n",
" if self.weights is None:\n",
" self.weights = np.log(np.ones(len(imp_samples)) / len(imp_samples))\n",
"\n",
" likelihoods = self.posterior.likelihood_function(imp_samples)\n",
"\n",
" unnormalized_weights = self.weights + likelihoods\n",
" self.weights = unnormalized_weights - logsumexp(unnormalized_weights)\n",
"\n",
" return self.weights"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we decide for a concrete type of model we want to test this on. We choose a Hidden Markov Modell (HMM) with discrete hidden states and discrete observed states:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"class HMMPosterior(Pdf):\n",
"\n",
" def __init__(self, trans_matrix, obs_probs):\n",
"\n",
" self.observations = []\n",
" self.obs_probs = obs_probs\n",
" self.prior = trans_matrix\n",
"\n",
" def likelihood_function(self, hidden_states):\n",
"\n",
" return np.log(self.obs_probs[hidden_states, self.observations[-1]])\n",
"\n",
" def add_observation(self, x):\n",
"\n",
" self.observations.append(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This guy takes a transition matrix for the hidden states (`trans_matrix`) and a matrix with the observation probabilities conditioned on hidden states (`obs_probs`). The likelihood function evaluates the likelihood of the most recent observation given a hidden state. `add_observation`, well, adds an observation. \n",
"The last ingredient now is a transition matrix object:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"class TransitionMatrix(object):\n",
"\n",
" def __init__(self, transition_matrix):\n",
"\n",
" self.transition_matrix = transition_matrix\n",
" self.hidden_states = np.arange(len(transition_matrix))\n",
"\n",
" def sample(self, old_states):\n",
"\n",
" return np.array([choice(self.hidden_states,\n",
" p=self.transition_matrix[hidden_state])\n",
" for hidden_state in old_states])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It basically just wraps a `numpy` matrix and provides a neat `sample` function that provides the required interface for use with `SequentialPrior`. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we have to stitch all of this together. First we define the parameters of the HMM:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"hidden_states = [0, 1]\n",
"observation_states = [0, 1]\n",
"start_prob = np.array([0.5, 0.5])\n",
"trans_matrix = np.array([[0.9, 0.1],\n",
" [0.1, 0.9]])\n",
"obs_probs = np.array([[0.9, 0.1],\n",
" [0.1, 0.9]])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's choose some number of importance samples which we will draw each time a new observation comes in:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"n_particles = 1000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now draw an initial state for the Markov chain in the HMM and instantiate all objects we need:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"start_state = choice(hidden_states, n_particles, p=start_prob)\n",
"\n",
"pdist = TransitionMatrix(trans_matrix)\n",
"prior = SequentialPrior(pdist, start_state)\n",
"posterior = HMMPosterior(prior, obs_probs)\n",
"sampler = SequentialImportanceSampler(posterior)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"No inference without some nice fake data generated from the model:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"n_observations = 30\n",
"# we have to start with some state...\n",
"true_states = [0]\n",
"observations = []\n",
"for _ in range(n_observations):\n",
" new_state = choice(hidden_states, p=trans_matrix[true_states[-1]])\n",
" new_observation = choice(\n",
" observation_states, p=obs_probs[new_state])\n",
" observations.append(new_observation)\n",
" true_states.append(new_state)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Slowly getting there. Here's two variables which will hold successively more samples and importance weights:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# will store sucessive posterior samples; is enlarged in each iteration\n",
"states_samples = None\n",
"# stores log-weights for all samples\n",
"all_log_ws = []"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Aaaaand we sample. We're in an on-line setting and pretend to have a new observation after each time step:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"for observation in observations:\n",
" # we're on-line: a new observation is arriving!\n",
" posterior.add_observation(observation)\n",
"\n",
" # draw samples from importance distribution and\n",
" # store them along with their corresponding importance\n",
" # weights\n",
" (samples, log_ws) = sampler.sample(n_particles)\n",
"\n",
" if states_samples is None:\n",
" states_samples = samples.reshape(1, n_particles)\n",
" else:\n",
" # append n_particles samples to states_samples\n",
" states_samples = np.vstack((states_samples, samples))\n",
"\n",
" all_log_ws.append(log_ws)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, to get samples from the target posterior distribution of hidden states, we have to draw unbiased samples again:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"# ... and back to actual probabilities\n",
"all_weights = np.exp(all_log_ws)\n",
"# get samples from target by reweighting importance samples\n",
"target_samples = np.array([choice(samples, n_samples, p=weight)\n",
" for weight, samples in zip(all_weights,\n",
" states_samples)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And finally, a fancy plot summarizing the results:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAn8AAADQCAYAAAB2mwfpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeXxb5ZX3v492yXvsJE5ik4UkhBAIS9gCaQk7DQVKKS2FaSm0TAtMZ8q8vF1oy7SdNm1Z+nYYOh1aGGgLZR8amhK2NCQhkIU4ISFxFm/xFnm3ta/P+8f1vci2ZMuObMvO8/187udKV8997pEsS0fnOed3hJQShUKhUCgUCsXxgWm8DVAoFAqFQqFQjB3K+VMoFAqFQqE4jlDOn0KhUCgUCsVxhHL+FAqFQqFQKI4jlPOnUCgUCoVCcRxhGY1JS0pK5Jw5c0ZjaoVCMUn54IMP2qSUU8fbjmxEfaYqFIrhMthn6qg4f3PmzGHHjh2jMbVCoZikCCHqxtuGbEV9pioUiuEy2GeqWvZVKBQKhUKhOI5Qzp9CoVAoFArFccSoLPsqFAqFIjlCiBzgN0AY2CClfHqcTVIoFMcZKvKnUCgUx4gQ4gkhRIsQYm+/41cKIQ4IIQ4LIb7Te/h64EUp5deAa8bcWIVCcdyjIn8KhUJx7DwJ/CfwB/2AEMIMPApcBjQA24UQa4AyYE/vsNjYmqlQKDJBbW0t4XCYhQsX9jnudrvp7OwkHA4TDAaRUhIIBLDZbOTm5mKxWPB4PAQCAePxeDxOPB7Hbrdjt9uxWDTXLCcnh9zcXEpKSnA6ndhsNoQQGbF/SOdPCPEEcDXQIqVckpGr9vJKRSMPvH6Apq4AMwud3HvFSVx3xqwxn0Mx9gz1dxur98ZY2JEJW8fquYzFHJMRKeVGIcScfofPAQ5LKasBhBDPAteiOYJlwC7U6otCkdUEg0EaGho4dOgQVVVVRCIRWltbWb16NQAulwu/35+RawkhsNlshvOXl5fHzJkzWbVqFZdffjmnn346ubm5GblWOpG/J+n3izYTvFLRyHdf3kMgov3wbewK8N2XtR/D6X5RZGIOxdgz1N9trN4bY2FHJmwdq+cyFnMcZ8wC6hPuNwDnAv8B/KcQYhXwaqqThRB3AHcAnHDCCaNopkJx/BEMBqmvryccDhOLxZBSAhCPx4lEIlRWVrJmzRp27dqF1+vF7XYnnWc4jp/VaqWgoACXy4XD4cBqtWI2m41onsvlorCwEKfTCUBpaSmzZ8/m7LPPZt68edjt9mN81h8zpPOX4hftMfPA6weMLwidQCTGA68fSPtLIhNzKMaeof5uY/XeGAs7MmHrWD2XdOcwESPHFCQQtxOIwIOv7+e6pdNBmNX/ZBpIKX3AV9IY9xjwGMCyZcvkaNulUEwG/H4/PT09AEgpDacONMdOSsnu3bt555136OnpQUqJEIKuri58Ph/xeByv18vBgwfp7u425hkOJpMJq9VKWVkZX/nKV7j99tspLS3N3JPMABnL+Rvur9SmrsCwjo/WHIqxZ6i/21i9N8bCjnTIhB367UJzD7mmAEJIgnEbTV1TALD7D7PAHsUkJOG4FXd0Cv64M63Xw9PTBrXP8BXn85xWfIhTHNXkmIN8o+47vNZ9IXOj78OzlwGwqVwQLxc82vJ5fuW+eci5JzmNQHnC/bLeY2kjhPg08On58+dn0i6FYkIhpaSnp8dw3hKP6/t4PE5zczOVlZWGQ9d/jJSSDz74gLfeeotwOIzD4SA/Px+LxYLf7ycajeJwOPD7/YTDYcrKyli4cCE5OTm89NJL2O12AoEAQghj3nnz5vGd73yHr33ta2P8qhwbGXP+hvsrdWahk8YkXwgzC51pXzMTcyjGnqH+bmP13hgLO9IhE3ZcPO0In3E8w1UFWzCLOABv9pzDv3VpeSkvLPgOxebOPue/3LmShzzf1+5svYOfzummyptHdyyP+fZ6tvlOYb3nHBZPCcCWW/mHYht7A/N4vvMymsJTqQzMBSBonwOn/hhkjKfercIXDLPdtzilrccR24EFQoi5aE7fF4AvDmcCKeWrwKvLli2bWN8sCkUaRKNRIpFIn2P9nTspJR0dHXR2ap9fJpOpj/Olj6uvr2fLli10d3cbhRL9HUCAqqoqTCYT5513HsuWLaOkpIRgMMiWLVtobGwkJyeHiooKAoEAHR0dfPjhh0ZEMBgMDpgvHo+Tn5+f+RdnlBm3at97rzipT24QgNNq5t4rThrTORRjz1B/t7F6bww6Rsb58UXwzrtv8HL7hXjjrhHZkQlb03kuPz3xBVw9Ffyu9TMcDpUTl4J2OZ17r9TGHJr9ID95v4pgVOI0hZhm7aA5VqbNEY+C+20+X9CIOT8EQDhuwR938F7wfL5w8cUwdxfragr5zv/uH2DHzZeuhFNvAaAo1sgvj8P/SSHEn4GLgBIhRANwv5TycSHE3cDrgBl4Qkr50TDnVZE/xYQjWXQuGe3t7cRiMcNJ6+/UgeYgVlZWEovFKCkpwWazGfMmjo9Go8Tjcc4991yWLFmC2WzG5/NRVVVlRAY3btxIdXU1HR0dtLS0sH79eiOq2NzcTFdXF6FQiK6uLuLxOEIIIpEIoVAo5XO54IILuO66647tBRsHRDpr2b05f39Nt9p32bJlMp0+lKra9/glG6t9Ty/28NPTK1hs2QNt70GkC4B/afkFfzl6CnMKzfzzFaeOe7VveaGVh887wLLAk7DiRcidC746/ro/wOq3mkb+mkrJ2g/28YcNH7CrPZ+SgvwxrfYVQnwgpVw28ldx8pLuZ6pCMd7EYjHC4fCQ4wKBAHV1ddhsNlwul3G8v5RJKBTC7XZTVFREaWnpgGVcgI6ODlavXk11dTVz586lsLCQeDxOR0eHUajR09PDvn37iMe1lRG9alYIQTweJxqNIqXE7/cjpcRqtTJlyhRsNhs9PT10dXUNcE6FENTX1zNrVnb6HIN9pg7p/CX+ogXc9P6iHewc9UGlmDDo738hoKMC1p0JBadAyXKYeoG2z5sP738FfHVw6d/Hz9ZYEA4/BvsfAv8RyD8Zzn0cpp4/fjZlEOX8pUZ9piomCj6fj+7ubiNCB8kjZl6vF6/XS0lJCTk5OX3G6k5WOBzm1Vdfpbq6msLCQsNJ7OzsNJaLm5qaqKiooK6uDrvdzrRp07DZbMTjcY4ePYrH4yEejxt7k8mE3W7HbDYPsCkYDBKNRrFarcycOdOoxO3s7KSjo8OwS3dQbTabsRScjQz2mZpOte9NmTdJocgSOnbA1q/C8qeh8DS4oQNsRQPHFZwMNU9BdyUULBp7OwEq/i8cfASmXgjL/hNmrQKhZOImM2rZV5EtJObmDRY0CgQCSClxOBx9xvfPv2tqaqKjo8MosEjGG2+8weuvv47ZbGbq1KmYzWZ6enqorq7G5/MRCoWoqakxIo0Wi4X29nbjfD3Kp0u56PMUFxcjhDA23e7Dhw9jMpkoLi7G5XIZzp6uu6ffN5vNmM1mTKaJ+/mrOnwojm8OPwaew+AqB5M5ueMHMO8r8OEP4PB/w1m/GlsbdXLnwsn3whm/HJ/rK8YcVfChyAaklESjUWDgsqyOEIJYLEZLSwuRSGTQThQ1NTW8++67CCE4cuQI8XicQCCA1+tFSonH42HLli189NFHRgWu7ox5PB58Pp+xvKzbZbPZcDgcfRw10LT12tvbsVqtzJo1i8LCwgH5glJKGhsbkVJis9kMR1OfQ59Tp6ioyHBEJyrK+VMcv0R6oO7PMPsLYCsYfKxjGpR/FqqfhKU/BYtr8PGjwaJvjf01FQrFcY++bKqLEieSGAXs6emhpaUFp9OJ1+tNWsBx8OBBnnzySbq6usjNzSUnJ8eQaenu7iYWi9HU1ER3dzdSSkwmE9FolEBAUzxIjOIlRvPmzp2b1OGsq6tDCEF+fj4zZ85MOqa7u5tQKITJZGLGjBlYrVbjuemt13RMJhMrV67kvffeG/kLmgUo509x/FL7Z4j6YP4d6Y1f8A2oexbqnoMTh9TozSydH0L+SWDOnMK7IvtRy76KsWCo3H+/3093d7fRlSLxHN3B6+npYfPmzbjdbmPZVI8YxmIx4vE4jY2NbNy4kYaGBhwOB6WlpYRCIbxeLw0NDYaock9PD7FYDLvdjsPh6NPT1mq1YrPZ8Hg8CCGYMmUKCxYsMJZgE5dydfFmIQQzZszoM0/iuPr6ekOyZcmSJYbt+uPxeBy3243ZbMZqtXLnnXcq50+hmLAcfkzL8ys+J73xU1fAub+H8jEu649H4O2VUHYtnPfE2F5bMa6oZV/FaKK3NhuKo0eP0tnZSVFRUZ8oWKLTuHfvXioqKhBCYLFY8Pl8eDwempqaiMVidHV1sX79evx+v+Fotbe395FTEUIYy7gOh4NFixaRk5NjaPvpWywWY/fu3djtds4//3xmzpzZp/pX3zc3NyOlpKioiMWLF2MymQbI0ESjUaMi2OFwUFtbO+D5t7a2ApoTWFRUxLJlE78uTTl/iuMTGdeWUa15WqVvOggBJ94+unYl4+h6CHdA2cTTklIoFNmLLl6sL6EmklhxW1NTg8lkMuRR9Md16urqeOyxx3C73ZSXlxs5efv27aOpqcmIDPp8Psxms5FDpxeRmEwmHA4HeXl5dHR0YLPZWLFiBXPmzMFkMhnFFUIITCYTNTU1fPTRR0yZMsUY07+AQwjBzp07kVJSWFjYp4+uEMJ4zg0NDYTDYUwmE7NmzSIajRqagbqjG4lEjMKP008/nbVr1xrL0BMV5fwpjk+ECebeMrJza56GQBMsvjezNqXiyPNgyYMZl4/N9RQKxYQmXZFlr9drOF7JKnJjsRg7d+5k//79FBQUGA5Pa2srsViMWCzG+vXr2bZtG+3t7YYmntVqJR6P09nZSSgUwmKxGDl1hYWFlJWV4XA4DEfMZDIZY9rb2yktLWXJkiV9cu4St4MHD+Lz+SgqKmL//v1Jn1skEuGjjz4yWsLV1NQYy8+Jr4le6JGXl0dRURHr169P+dqFQiHWrl3L3/72N6SUFBQMkSuexSjnT5HdBNzgnJ7ZOSMerWp33lfAXjz8891vaw7Zgn8E6yi39YlHoOF/tSVfs2Po8YpJhcr5UwyXZC3TkuHxeGhpaSEnJydlZW5bWxvbtm2jo6ODkpISfD4fra2t7N27F4/HQyQSoaKigmAwaDiPfr/fWF4VQuByucjJyUFKSUlJCXfccQcnnHBC0utt3LiRuro6zjnnHAoKCjCZTEm3TZs2YbPZuOqqq/oUeiTq8DU1NbF9+3ZKSkq45pprDGdTH6fzzjvv0NLSwt13301XVxdvv/02S5cupbi4GIvFgtlspqamhgMHDjBz5kwuvPBCLBYLJpOJ66+/frh/nqxBOX+K7KX2GdhyM1yxHYozmGNR92eouFfTyxuJ8zf/61D9P1D7tFYEMpq4N0C4E064cXSvo8hKVM6fYriEw2F8Ph82m21A5CqxSMPr9QKQn59vaPL1lz959NFHqaiowOFwMGvWLGw2G9XV1Rw5coRIJEI4HDbmcTqdFBQU4HR+3MdbCEFRUREdHR3GsmplZSWVlZUD7JZS8s477xj6fboEjB6t07empiYOHDhgRCU//PDDpK/DgQMH8Pv9mM1mXnvttZSvV2trK2azmTPPPNOI+v3qV79i5cqVxphf/OIX/PCHP+TCCy/k2WefHfwPMEEYf+evcS1EumHOsPqdK44HGl/V9t17M+v8HX4MCk+F4nNHdn7x2VB0Bhz6L80RTDdncCSUXgqXv6ddT6FQKIYgFAoRDoex2WxJRYiFEHR3d7N37158Ph89PT2YzWaCwSA+nw8pJd3d3fz9739n48aNRKNRSkpKcDgcmEwmurq0tpcul4upU6cSDoeZNWsWd911F2eddZaRf6dfC2D16tU0NDRw5513MmvWLEPCRS/S0PX7tm3bRnFxMaeddpphv57vp28VFRXs2LGDsrIy/vVf/3XAtfT5fv7zn+Pz+fjWt77F2WefbcyjCzTr9//0pz/xxhtvcO211/ZZ8p3sjL/zV/U4eA4o508xkFio7z4TdHygbcv+c+ROmxBaxG/bHVoP4KnLM2dfsmuVnDd68ysUigmB3ns2kWT3Ozo6iEajfSJw/XMAN27cyJ49ezCbzbS2tiKlpL6+nnA4jJSSDz74gP379+P1eo38PY/HA2h5gpFIBKvVSmtrK9FoFJvNxr59+5Lm30UiETZu3IjJZOLZZ5/t0/kjkY6ODmpqasjNzeWDDz4w8v30yJ9elfzBBx/Q09NDc3MzP//5z1O+XjU1NTgcDj71qU9x9tlnpxy3ceNGnE5n0nZvk5nxd/5c5XD0rfG2QpGNlF0LrjItty5THP4dmJ0w5+Zjm2f2Tdqyb3zoBuYj5uh6qH8RTvvJyJanFQrFpCAejxt5fIN1zgiFQvj9fqSUBIPBPjlwugPY0dHBjh076OnpYe7cuYYzZrVasdvtlJSUcOjQIWw2G+eccw4XXHAB119/vVHc8Oabb3L06FFyc3PZvHkzAA888MCAvregOZotLS2sX7+eJUuW8KUvfalPBC5x7JYtW9i6dSsXXnghl1xyyYBInX770KFDRKNR/umf/olFixYZ+XeJET2z2cyDDz6IxWLhrLPOyuwfY5Iw/s5fTjlEPRDuHrrLguL4Yt6XtS2ThNo0x81WeGzzWHPh0g0ZMSklNX+AhlfgzP83utdRZC2q4EOhE4lEUi7l6o5dZ2cnR48e7VOZ25/nn3+eLVu2GPIrLpcLr9fLgQMHaG9vx2Qy0dbWZuQO1tbW8pe//KWPjp7dbsdms3HkyBHy8/N56KGHUtq9d+9e3G43drud559/3oji6Z06dBoaGjCZTJxzzjmcd955Ro6fPk6/bbfbKSws5Jvf/CaFhak/x59++mmACd1/dzQZf+fPVa7t/fXK+VN8TDwG0R6oex6CLXDqDzIz74oXtbkzRbgLvDUwJcM5ebGQ5viVfwbMtszOrZgwqIIPBWhFHB6PB5vNlnJ5MhaLcejQIRobG7FarUYnDl2rTt93d3eTk5PDNddcw5VXXkleXh67d+/m+eefp7W1lXnz5vHqq6/S2dnJxRdfTGlpKRaLxah8tVgslJSUcOTIEV555RVWrlzJ1Vdf3Sc6l7j/9a9/TUNDA5deeimnnXbagJw7s9lMOBzmvvvuIxKJ8NJLLw1aoNHS0kJhYWHK5WNFemSX81e4ZHxtUWQPvhp4dYF2O2duZpw/XTbGlMHcjk03gP8IXF2paQdmiqNvaoVQqspXoTjuCQQCeDweioqKsFgsSXP9Dh48yDvvvIPf76ejowOLxUI4HDZapvn9ft5++20aGhqIRCK0tLTwwgsvAFo7Nb3bRmNjIz6fj6lTp3L33Xdz0kkn9am21aNwDz74IMFgkHA4zL59+1Jq8jU2NlJcXMznP/95li9fnnSc2+0mNzcXp9PJz372M6ZPn55S5uVrX/sa5557rnL+jpHxd/6mnAXXNYIjw1puiomN57C2n3oBtG/TOnIci3PVsRNePxtW/C+UXZMZGwHm3Qrv/QO4/w6ll2Ru3rrnwVYE0zM45yRGaMlDNwPzpJQ/FkKcAJRKKbeNs2kKRUoSl1JTPQ7acq7P58PpdBoRvMR8OSEEoZBWGHfeeedx/vnnY7FYaGxspKGhAbvdzoEDB+jp6cFkMmGz2XA4HESjUcOpslqtWCwW/H6/0WrtmWeewW5P3k983bp1hMNhotGoIeCsizXrtwGampqIRqPcf//9STuJgBbZbGpqYs6cOZx++ulMnTo15Wvicrn6FLMoRsb4O39mO7hmjrcVimzDW6Xtp10Ere9CsPXYxJ4P/w5MNpi2IiPmGZxwA+z8F032JZPOn2Oa1kpOLfmmy2+AOHAx8GPAA7wEpC7zUyjGmVAoNKjzB9pyrtvtNhw1vbCify/bmpoaOjo6CIVCHD16lHg8Tm1tLQ0NDSxYsIBQKEQ0GmXhwoWcdNJJ3HfffUYe3YEDB9i+fTt2u53u7m7eeecd7rvvPk444YQ+bdX0TQjBvn376O7u5o9//KPRrq0/gUCA3/72t8ydO5ebbrrJOL+/fIvX6+XBBx/kc5/73KCOnyJzjL/zB3DgEa3H6rxbx9sSRbbgOQxml6apB1pawEidv4hXq8w94UYtmpZJzA6tU0jlr8DflLkfMmc+mJl5jh/OlVKeKYSoAJBSdgohUnrOQojlwBwSPgOllH8YdSsVigT0YgbdeUpWLRsMBvH7/dhsNkPupL+2XUVFBc899xzt7e00NzdTUFBANBpl3759HD16FLPZTE9Pj9HVY//+/bz77rsDHM+pU6fS2NhIJBJh3bp15OXlDZBckVLi9/uprq7G4XDwwgsv9BFiTrytLwuXl5dz2223pXwd2traePLJJzn11FMz+voqUpMdzl/tM2BxKedP8THeKsg7UcsJteZrXS5GSt2zWkX5/AxKxiQy/x9h/4Nay7dF/3Ls83lrIGfO6IpHTz4iQggzIAGEEFPRIoEDEEL8ETgR2AXo1T8SyDrnT1X7Tm66u7sxm81GcUYyAoEA0WiUsrIyZs+ePUDqxe1209bWRiQSYf78+Vx00UUUFhbi9XqN5dgFCxawZ88ejhw5wsKFC1m+fDkOh6NP8YXT6SQSifDMM89QWFjIqaeeis1m6xP10yN2O3bsMHr2Pvfcc32ieImFHLozW1ZWNtovpWKYZIfzl1MOnbtGfv7R9VCwGJylmbNJMb7M/RJEfVpni891H9tcVY9DwSlQcn5mbOtP3nz4XI8WvT5WYkH421KYf4eK/g2P/wD+F5gmhPgpcAPw/RRjlwGL5VDrbVmAqvadvCTKpuTm5gIDNfzC4TA7duzgww8/xO12c+DAAaSUtLe3E4lEiMfjvPrqqzQ0NNDc3MyMGTMMSZhgMEhNTQ0+n8+QbsnPz+eqq67i5ptvThqpO3LkCLm5uZxzzjnceeedKW232+2sWbOGr3/96326bPQnEAhw3XXXcf75o/TZqxgx2eH8ucq1Vl5SDj/aEfXB+ktg2ifg0ndGxz7F2HPCDZmba8WL4G8c3UhaJhw/gObXtSjljMszM99xgpTyaSHEB8AlgACuk1IObDegsRcoBZrHyj7F8UOiIHN/En9v+P1+Wltbyc3NNTpr9B/T0tJiCDI7nU6EEITDYWpra43cO72QIxwOEwgEcLvdWCwWotEofr+fWCxGW1sbra2tOBwO3n//fTo6OpLa19bWhtvtxu12s3bt2j6afImO4vvvv08kEkEIMajotCJ7yR7nLxaEcMfwOxl4a7V9y8aMm6UYJ6I+8FRB/kItp27nv2pLv6feP7L5XLO0bTSpfUarKD7WaF3d89r/wPSVQ49VGAgh/gN4Vkr5aBrDS4B9QohtgNE7UEqZwTJwxfGK7iDp1a79q3J1dAfR4XBgtVqTjnO5XNhsNpYuXco111xDTk4ONTU1vPXWWzgcDkP3z+l0UlRUxNKlS8nPzwcwIn0mk4nCwkICgQCf+tSn+NrXvmaIRfdfqv2f//kf3n//fTweD3/6058GPDd9ibi+vh6AefPmZfjVU4wV2eP8CQsEjg7f+fPVfHw74tU6LygmNu3b4e2VcPGbUHqp5lTFwyNz/nSR6JPuzrydibRvh8OPwRkPjDzCGA1A4xqtA4kpdQ6QIikfAN8XQpyEtvz7rJRyR4qx/zZmVimOS4QQKSVSdGw2G3a7neLi4pSadevXr2ffvn00NTURCAQwmUxs2LCB3bt3EwwG8fl8NDU1GdcLBAJGnp5ux5QpU4jFYjgcDi644AJOPvnklDa5XC4sFgs/+clPOPXUU/u0TUt0Tl966SX++Mc/smrVqhG8OopsIDucv7Jr4fPBkYnvehOcP38DFCzKnF2K8UHX+MvtTXJ3lUPLCJf0a/6gvS9G2/lzlUPMrxWm2KeMbI7mdRD1wmwl7DxcpJRPAU8JIaYAnwV+IYQ4QUq5IMnYd4QQ0/lYBmablLJlDM1VTGL0/rq67h70XcrVb+viy8nw+XwcPnyYxx57jIaGBnJzc6mqqsJsNlNbW0tHRwcOh4OioiI6OjooKCjgsssu4xvf+AaFhYVYLBZDt89qtfLII4+wdu1a3n33XSoqKvos5yYu61ZWViKEICcnR2npTXKyw/kzHYMZfZy/euX8TQa8VVrkS+/+klMOgUatLdtwfyB4q6Eg9S/djJGjd6ppGLnzN+NK+MQrmrahYqTMBxYBs4GkOX9CiBuBB4ANaPmBjwgh7pVSvjhWRiomL+FwmFgsNqANW/9lXZtNUyLS94m88MIL/PnPf2b37t1IKYnFYkaOX09PD4FAwHDgXC4Xn/zkJ/n973+fUm/v0KFDxGIxpkyZQlFRUZ/2aroos34sFouhqssnP9nh/AFsvxuKlw1f7mXhnZC/ALbfqX3xKiY+nsOa1Inu6LnKQMYgeHR4uXsyrjl/s8ZgaSKxTWHRaSObw+LUouCKYSOE+CXwGaAKeA74iZSyK8Xw+4Cz9WhfryzMW4By/hRJkVISiUSSCjL3PxaNRjGbzUYFbyoaGhrweDyYTKYBgs1bt26lrq7O6LxRWlrKnDlzsNlsWK1WWlpaKCoqIhKJYDabKSgo4IknnhhQnKFv1dXVlJSU8NWvfpXS0tSqGCaTiY6OjiGXrBUTn+xx/prWQqRr+M5f3nzImQ35i6BACUROCryHP17yBe12wRKIeIY3T6AZ4iHIPTGz9iXDVQaWvOHbqBNwQ9XvtHy/vDGwd/JRBZwvpWxLY6yp3zJvO5DBxsyZQ+n8ZQeJ0TdIXcQhhMBisaRVAevz+YxcPSEE9fX1dHV1EYvFaG9vJxwOM2vWLAoLC7n88stZsmQJsVgMp9NJR0cHJ554Ips3b8Zms1FaWkpLS0sf3T6z2YzdbsdsNpOTk8P8+fMpKSnJ/IujmJBkj/PnKteiJsNBSqh5StNvU9WRk4czHgCRUPAw4zJYtWf48/jrtX7AuWNQkeacBTf2jPz87o/gwx9ovYyV85c2QohFUspKYDtwQm9PXwMp5c4kp60TQrwO/Ln3/ueBv42upSND6fxlF1ar1XAABxuTroSky+Vi+nStc1FFRe5d1tMAACAASURBVAVr1qwBYOfOnbS2tuL1enG73QSDQTZs2ABAY2MjdrudWCxGd3c3J598MvffP3gx3MGDBznzzDNTLgsrjj+y553gKoe294Z3TrgD3v8KnPmwtuQb6Yby60fHPsXYUXppZuYpOQ9uDIxNp4xjvYavVtvnzDlWS4437gHuAB5K8phE6/Xb96CU9wohPgtc0HvoMSnl/46eiYrjBSklXq+XeDw+4Hj/+5FIxCiqaGtr4/3336e5uZny8nLC4bDRSk0XYNY19YqLiykqKmLu3LlUVVUxe/bsMXt+islD9jh/OeVQ/4KWpyXSXIHRiz1y5sLBRzRtOOX8TWz8TVq3l2kr+gon//0q7dgp3xvefOaU7V0zz54faz9Izvp/wz/XV6u9712qDdJwkFLe0XvzKillMPExIURy/QztvJeAl0bTNsXxRygUwuPxIKUcECEUQhjOnJQSn8+Hw+Ggs7OT7du3c/jwYTo6OigrK8NkMpGfn8+0adOYMmUKV1xxBfPnzzeKPk477TRmz55NR0cHy5cvH6dnq5jIZJHzNxecMyHclX61pK7xlzsXnGXg3jBq5inGCPd6eO8fYNX+vpXbvhroHGYXjX0PaOLhp/4gszamonsfdOwYmfPnrdXew0rfb6RsAc4c7JgQYrOU8kIhhIfeHsD6Q4CUUuaPvpmKyUYyGZeSkhKKi/tq1h45coTHH3+cWCxGKBSivr4em82Gw+HA7XZz8OBBzGYzVVVVhtxKYWEh5eXlXHHFFZxxxhl95otGo6P/5BSTluxx/hb8o7YNB2+C85dTri37RjyZa7WlGHs8hwGh/U0TcZWDb5g5ofUvakUYjJHzl1MODa+MrE2hvx5y54yKWZMZIUQpMAtwCiHOQHPkAPIBV+JYKeWFvXv1AaEYFpFIhM7OTqMzRiqCwWDK1m6dnZ3EYjEuuOAC3G438Xic0tJS8vPzqa6uNuY/6aSTaG5uRghBPB6noKCAsjK1IqDILNnj/I0Ebw3Ypmitv5y9/xz+hrHRdVOMDt4qzdEz95MacJVB85vDnKt6bNMAXOVadXGoDRxTh3fuxW9pPX0Vw+UK4FagDHg44bgHSJojIIT4o5TyH4Y6plDo6Eu1DodjUP0+k8lEMBgkJycn5VznnHMOb775Jrm5ucyYMYPi4mJ8Ph95eXlG5W5OTo7hRE6ZMoWCgoLReWKK45bscf6ifth0A8y9BeZ8Mb1zlv4UFt6l3U4U2VXO38TFc1iT7+mPqxyCzRCPpicKHunRnLCxkHnRSdT6G67zZzKDrTDzNk1yEjp7fLY3jy8dTkm8I4SwAGdl3DhF1pOor5d4rD+xWAzQ+vAOVjHr9/vp6ekxumskzud2u+np6aGpqYm2tjZmzJjBpZdeSmlpKQUFBRw6dAi73c5ZZ53F1q1bCYfD2O12Fi1alFQIWqE4FrLH+TM7oeXvULA4fefPPuXj/MApy+Dqg5rmn2Li4q2CsusGHi86HUqv0NqfpeMkeaq0/VjIvOjkztXev7HQ0GMT8TfCnh9pP2SKlo6ObZMcKeVLQohVaI6dI+H4j/XbQojvokUDnUIIXZdHAGHgsTE0V5ElhEKhtGRZ9Py6ofT79O4bkUikT39dfS+EwGQyYTKZcDgcg0YIFYrRJHucPyG0pdt0tf5kHPb+BGZdDVPOAotL6/ShmNhc8jaYkhRpll8/vCXcSLdWQDSWkb+i02HVR8M/r+dAr8DzFzJv03GCEOK3aDl+K4HfAzcA2xLHSClXA6uFEKullN8deysV2YaUErPZPKT+nRCCWCyWlnizxWKhvLwcq7Vv8VZbWxt5eXlMnToVu92OzWYjEAgQjUbp6ekhGAwSj8fp7u4mFAoRDodVxE8xamSP8wfa0m26LdoCTbDn38AxXXP+AA7/DqwFMPvGUTNRMcoUZqhLy/SL4DONmZlrtNE1/lTBx7GwXEp5mhDiQynlj4QQDwGvpRi7TQhRIKXsBhBCFAIXSSlfGTNrFVmDHo0bDL/fj8fjMcYmtmLT0TX+hmLXrl28++67WCwW3G43JpOJ6upqDhw4gNlsRkpJTU0N8XicxYsXK2FmxaiQXS2NhtPlI1HjT+fwf0P1kxk3SzFGdO6GQ/8NkSQfoFE//GUOVI5ARmUsefdm+OCe4Z3jq9M0/pyqou8YCPTu/UKImUAEmJFi7P264wfQ2wN48BYJiuOa7u5u/H4/oVCIUChEJBIhEokQjUaN/rlSSoQQQ+YFbtq0ierqajo6OnC73bS2ttLT02PMm7gUPW/ePE455ZSUcykUIyW7flIULgXfkfSkMhJlXnRcZVqFp2Ji0vQ32P295DmfFheEOj7+uw/Ftm9oVeBn/CKzNg5FoAH8R4Z3jq9Waw83loLUk4+/9kbwHgB2oun4/T7F2GQ/esfss1AIMQ+4DyiQUt4wVtdVHBs5OTlMnz590ChhNBrtk++XjFgsRm5uLtdddx2XX345drudbdu2YTKZsNvtfOpTn6Krq4t4PM7tt9/OokWLUs6lUIyU7Ir8nXwPXPr39DTSfDWA6FvgMRItOEX24K3SlvFT6TTmDCMyfPSN4TthmcA1jNQFnVgweYWzIm2klD+RUnb1VvzOBhZJKVMJPO4QQjwshDixd3sY+CCd6wghnhBCtAgh9vY7fqUQ4oAQ4rAQ4jtD2Fotpbw9nespRpd4PE4sFiMajRqbHtWLRCKEw2HC4bBR7ZspTCYTeXl5FBYWUlBQQE5ODjabDZvNhtPpxGKxYDabVc6fYtTIrsjfcPAd0RL6E/XgXGUQ6dKWDa2542ebYmR4Dg9eoJFuQVA8oi2ljkcBhascAo3Da1N44XNatFsxbIQQKauAettpvZzkoX9CU/5+rvf+m8BdaV7ySeA/gT8kXMcMPApcBjQA24UQawAzsLrf+bdJKVvSvJZiFJFS0tXVhcViGVCc0Z9QKITNZhs0otfU1MS6desIBAIUFRUZ2oD6knBdXR0VFRWGw6nIXvQuLOFwmGg0it/vH2+TMk52OX++OthwNSz9dyi7dvCx5/4Owp19j+k6a4EmsC4cHRsVo4e3CqZfnPrxnHLo2j30PP56kLGxrfTVcZVrzmfQDc5UKWdJGG5HEIXOpwd5TAIDnD8ppQ8YNDqXckIpNwoh5vQ7fA5wWEpZDSCEeBa4tre6+OqRXKd3njuAOwBOOOGEkU6jSIGeV2e323G5XAOkWRJv+3w+CgoKBnX+GhoacLvd5Ofn43Q6EUJgNpsNaZdgMNhHrHnp0qUqsjeOtLe3c8MNN1BZWUlnZycLFy7k6NGjAHzhC1/Abu/baCCdSu+JRHY5f9YC6N7b2+JrCIQJ7H17J1J+PdzoB4tzdOxTjB6xoLZcmjvI8uf0i0FYhs4JNTT+xsH5KzgFZlyRvtafvwG23gFLvg9TVYP24SKl/MpwzxFCLAT+DzCHhM9AKeUgvzwGZRaQGJJuAM4d5PrFwE+BM4QQ3+11EgcgpXyMXv3BZcuWqdBwmvSPrKXS8dPHWCyWIZ0ws9mc1pe/zWbjwgsvZOnSgXqdu3fvprGx0YgIzpkzZ8j5FKNHfX09Bw8epLS0lFgsxs0338xrr73G1q1bufHGG1m8eLEh1m21Wlm9Oum/6YQl+5w/S+7QS3uxMGz/Bsz7Mkz7xMfHzUn04RQTA7MDbuhEC9akYM4X0xMAFyYoPnd88uimf1Lb0sVzGJpfg5P/dfRsOg4QQvww2fFEkecEXgB+i1YQMubrb1LKduDr6YwVQnwa+PT8+SonNF30Ctxkzlr/YyaTaUC7tpEgpTScuuGiO6sjOVdx7CxbtoxoNMr9999Pa2srW7du5XOf+xwrV67sM+6RRx4ZJwtHh+xy/oRIT+7FVwfVT8C0FX2dPymh4v9AyXI44bOja6si86TTuUPGtSVd0yA5OqWXaNtEQNf4y5kznlZMBnwJtx1oy637U4yNSin/K4PXbgTKE+6X9R47ZqSUrwKvLlu27GuZmO94QZdcGYx4PI7dbh9S4y8cDtPU1ITL5SInJ6dPS7jE201NTYTD4T4OZnt7u6HZd/DgQRoaGojH48TjcTZt2mScu3fvXqLR6KRbWlRkL9nl/IFWtDFUxa5Pl3np17pLCKh9GiIe5fxNNJpeg9Z34dT7Uzt2nsPw15Ph/KfSbwE4HqxdArM+DaensUzgrQXEx/mqihEhpXwo8b4Q4kHg9RTDXxVC3An8L2Csz0spO0Z4+e3AAiHEXDSn7wtAFr9BFcNBT/wvLCw0cvkSN9CczYKCAnJzcyku/jgdqauri0AgwNSpU8nPzycej+P1eonH47S0tBjOY1tbm6H1t23bNlpaWvrMoxge4XCYyspKQqEQ999/f5/qbX3vdrvp7Oxky5YtQ/4AmIxkn/NXeolWyTsYyQSedVxlw5faUIw/jWuh9k9w2k9Sj3GUgowOHRl++2IoOBWW/TqzNqZLPJy+3qSvFlxK428UcKFF4JLx5d79vQnHJDBkI2ghxJ+Bi4ASIUQDmmD040KIu9GcTTPwhJRyBH3+kl5PLfuOA4mRPT03cMqUKUyZMiXlOXl5ebhcrgGFAiaTiUWLFhEKhQgEAkZ0MBKJGMvOVqvVyDuUUmIymSgoKGDatGmj9AwnNx999BEfffQRZrOZtWvXGjI6VqvVeK2DwSAAZWVlXHHFFeNs8diTfc7f4m8PPcZXo0WHnDMHPuYqT69gRJFdeA9rOXqDLXtYc8FaOHhkWEro+AAKlmTexnQZTqcaa4GWn6g4JoQQe/g4YdQMTAWS5fshpUzyqzE9pJQ3pTj+N+BvI513kOupZd9Rwufz0draitfrxWazJW3ZBhgt2zKxJBsIBDCbzZSWlnLGGWcY3UE8Hg/19fXE43Fyc3PJzc1l3rx55OYqybJEotEoXV1d1NbW8tZbb/WJ5CVuhw8fJhqNsmLFCtauXZt0rl27drFq1Sq+/e1vc/HFI631mrhkn/OnM1hFZ9QPeQvAlCRR11UG7g2jappiFPAc/rhH82AMJfQcaodIz8CUgLHEVQ7ut9MbO17RyclHoqRKFHBLKaPJBgohvpTsuJTyD8mOjycq8jd6hEIhI+9Pzw/sv5Sry7VEIhFcLtcxXa+7u5v6+nrMZjMej6ePU1JTU0NNTY1xTb/fz9SpU4/pepORNWvWsHv3bg4ePMi2bdtSjgsGg8TjcVVRPQjZ5/y1boENn4JPrulbzJHIsv9ILYrrKgeTBaIBJfkyUTBEmT8/9Nihomr6cut4yLzouMo1rcl4VHsvKkYdKWWdEKIIrfDCAkzvFXnemWT42Qm3HcAlaC3hss75U5G/4ROLxQiHw8b9/hE9fe/3+zGbzRQWFg5Yqk0kHo/jcDiOOS8sEokAMG3aNBYuXMiXv/xlzGYzZrOZHTt2sHHjRkwmE8uXL+dvf/sb1147hNbtcYgehb322mu55557sNls2O12YylX3yoqKrjpppv44hdV6m0qsu+byVYEke6h8/ZSRQVPvhcW/9/M26UYPYJuTeolHYdt9hch1Jr6ca+u8TeOkb+py2HeVyAWAFOKVnWg5bZuWAVnPAgzj7+ck0wihPgJcCtQxcfLvxIYsJ4jpfynfucWAs+OsomKMcLn8yVtydZfxDkej2OxWLBYjv1r0Ov1sm/fPg4fPszmzZtxOp1Eo1Hq6upoa2ujtraWiooK2tvbKSkpGVAYYrfbM2LHRGX37t1UVVXxm9/8hpdfftlYvtWLbfR9c3Mz0WiUGTNmDNrzWFVND032vdv0qsdU0Z2IBzZ9FhbdAzOvHPi4+qNPPFxl8LkeTcJlKObePPjjjulQ/lnIHXFa17Ez8yptGwpvjSZqnix9QTFcbgROlFKGhxw5EB8wjm8YRTok6ugl5uX1v607dUVFRUBqR0DX18uEzt+uXbt48cUX8Xg8uN1uo2OIx+PB5/PR3t5OQ0MDwWAQk8mkOrb0Y/PmzfT09BCPx8nJyWHKlClGFM9ut/eJ6L3xxht84Qvj0LpzkpF9zt9QSf3eGjj6JpyYoi96uFPrmDDvVpi1atTMVGQYIbTuHUMh41qk0FaUXNS79GJtG29kXFv2HayKV2n8ZZK9QCEwZN9cIcSrfBwdNAGLgedHz7SRo3L+NPov5SYj0cmzWCyjHv3RncdYLIbb7SYej3P66adz1VVXUVJSgsViobm5mc7OTs4991zef/999u3bx+WXX87tt6f4/ppkbN26ldraWl5++WX27NljRPESI3qhUIi9e/diNpu5++67ufrq1B0R//jHP7J9+3bKylIV8ivSJfucPxg8qd83iMwLgNkF9S9C4anK+ZsoHHwUuvbCOWno7ro3wPpL4JL1MH3lwMdjwfHv9BLxwktTYOnP4OT/k3qcr07bK42/TLAaqBBC7KWvdt81+m0hhF1KGQIeTDgvCtRJKbNSH0rl/PXFZrP1yb1L5uD1F1o+Fnw+H++++y6tra0cOHAAIQSxWGxAN49du3YRCASYM2cOy5Yto6SkBNDkX44ePcqJJ55IbW0tNpuN/Px8rNZBROonAK2trTQ3N/PCCy8Qi8UGOHP6/VdffZW2tjaqq6sRQhhRPLvdTk5ODna7HbvdjtfrJRqNsny5anE5VmSn8zfnltRf4LrGX6plPbNdW/pTWn8Th+Y30tfFc/X+4ksVGX51Acy6Bs5+NOnDkUjEWH4ZVea/CuFc2J+qyQQQuxhOOh8OpvncJwkOh4OysrJMfwE+BfwC2APEU4x5DzgT+KqU8h8yeXHF2JBYjTsWdHV10dXVhcPhoLi4GKfTaWjzJW6dnZ1UVlYOmoc2mXjzzTepqanhySefxGQy9XHqEm/rGodPP/30oJW3Dz74IMFgcFAdRUVmyU7nb7CCDW8NWHLAXpJ6jKssfZ01xfjjrUq/D6/u/CX7+0YDmtPvmJ7y9IaGBvLy8pgzZ87ofol0xbQfMIM9r0CzJgidM3v07MgypJRG/tPcuRlNs/NLKf9jiDE2IcQXgeVCiOuT2PZyJg1STB4WLlzI4sWLsVqtRl5h4qZH85zO40NhIhqNkpeXx8svv0xubm7Kz9I77riDoqKirJVckVKyadMm6uvrCQQCVFRUAPC73/2OV155hUAggN/vx+fzUVtby4knjqOKRIbJTucPtKUzs3NgMrw1D6auGLywQwk9TxxkXHP+ZiQp3kmGxQX24uTOn55DN0jVcDAYHH3HD8Bk0yRsBsM5Y3RtyEKEEBQXF9PaOkjF9sjYJIRYDayh77JvotTL14Gb0XIDP93vfAlknfOncv7GH11UGEhZHHL06FF8Pt+gYyYb+jLuRK6s9fl8/OxnPzP+Znpe6dq1a8nPzzf6OTudTk4//XRuuOGG8TQ3o2Sn81f7DGy5Ga4+CPkL+j629N+HPr/gFAh3jY5tiswSaNby9PKG8YsqldZfmjIvY/JhZbJpkkWpGEzEfJIzSq//Gb378xKO9ZF6kVJuBjYLIXZIKR8fDSMyjcr5G390h6CoqMhYyuy/TZ8+nby8PObOnWtUGU80pJTs3LmTqqoqfvvb3xIMBgmFQgSDQWPT7x84cACY+JIqet5mZWUlLpeL++67j8cee4xXXnmFlSuT5JRPIrLT+dMjIv76gc5fOqTjICqyg3An5J8MeSelf87J39YigP3R8waH40iOFrZCzQFM5eTFw9C9T1vytas8l2NFSpn2J/VEcfwU2YPFYmHGjBkp81SdTicWi2VQsehsp7m5mfXr19Pa2sqWLVtwOBw4HA7sdjtOp5Pc3FzjWFNTEzk5ORO+cAU0B1aX3pkMzyddstT5S5HXFWyDt1bA0tVQft3Y26XIPIVL4Op9wztnTgqNp6IztOpaexa0RbIValsq4mFN1zAdeRvFkAghfpjsuJQyaX9fhcLj8RiFX4l5fIm3m5qa8Hg8x3SdeDxOKBQiEAgQi8Xo7u6moaHB6EsbiURoaRlSoWjExGIxampqCAaDhMPhPpG8YDBIIBAgGAzS3t7O0aNHOfXUU/nDHwZvdlNXVzdARFsxscjOb55USf3eauipBDFEm53uffD+bXDGAzBtxejYqBg/Ij3QvV9z9hJ19KatyJ6/t5Qgo4ApuYhzvDctbTAdQMVw8CXcdqD1+h2k1HpioHL+RgcpJW1tbUZvX5PJhBDCWMY1m82YTCZyc3NxOp0UFRUNOyrk9Xp54YUXaG1tZc+ePVRWVtLZ2cmOHTswm83EYjGi0aixSSmNHsOZpK6ujsrKStra2pBS4nQ6jQhebm4uU6dOxeFw4PP5qKqq4rrrVGDleCA7nT+LU6vm7e/8+YaQedExO6B9q1b0kS3OgCI5O76p5cad/1T65zS8Cu/dAldXQn7CcrG3Gpwzx1/nD7Q8xu6PtPeqvTjJ472Ctaa+zl8gEODKK69k/fr1GU0cv+222/jrX//KtGnT2Lt3r3H8oosu4sknn8xoNV44HObSSy9l/fr1Y9aySkr5UOJ9IcSDwOupxgshZgGzSfgMlFJuHDUDR4jK+cscUkojqqdr9Q1ViRqJRMjNzSU3N3fY1/N4PIRCITo7O6mtraWpqcm439ra2kcyxmKxGA7nsTqAUkrC4TCBQIBAIEBdXR1SSq688kp+8IMfpDyvoaGBjRs3csoppxzT9RUTg+x0/gBOuW9g1aZ3CIFnHecsba/kXrKf1nfBMW145yRGhnXnT8Zh7Smw4C4488HU544VulMXT9GVIB4Gk3VAFPuJJ57g+uuvz3jF4K233srdd9/Nl770pYzOmwybzcYll1zCc889x803D9GOb/RwAUnbAAghfgF8HtgH6GtXEsg6508xPGKxGB0dHYBW2a8v4+qbTjweJxKJjEnBQiQSYeHChZSUlFBRUcFZZ53F5z//+T56ePre6XRy6NChtObdv38/u3fvxu/3s3nzZoLBIH6/n2Aw2Oe59vT0EI1GJ3Q+oiLzZK/zt+hfBh7zVmsRQesQv8LMds2hUELP2Y2U4D0MJecP77yc3o4YiULPetXwEJW+A3jrooHHTrgRFt4JUT9s+NTAx+fdqm3BNtjcr/T/0g3a3mQGYU4p93LTbfcQj0epOXIUt9vNb37zG1atWsXTTz/NM888Y4xbuXIl3/ve97jsssv4/ve/T3d3N4888sjwniPwiU98wpCrSMXevXu544472LJlCwA7d+7k3nvv5e233x729a677jq++93vjpnzJ4TYw8ct28zAVCBVvt91wEm93T4Ukwi9y4TdbjciaolLuvpeSonH4xmzytx58+Yxd+5ccnNzmTt3LmecccbQJw3B9u3b6ejowGKxUFpaisPhwOl09tkcDgetra0Eg8FJX72qGB7Z6/xF/eA7AgUJiul5C6DsM+mdn0oORJE9hNq1/L10BZ51kkV29UrfQTT+xhyzLWXkb/feSq699lqee3ENmzdv5p577uGyyy6jurq6zzLUj370I374wx/S0tJCRUUFa9as6TPPihUrkiakP/jgg1x66aXDMnfx4sVUV1cbze7vueceHn744bTP7+zsNL5MlyxZwvbt24d1/WMksSFoFHBLKaMpxlYDVhL0ABXZj54jJ4Tos4SbGNnz+/3EYjEKCgoGdezi8fiAVnETEavVyg033MAVV1yRcszhw4fJy8s7ripZFUOTvc7fwUdg13fgcz2asDPA4nvTP3/aJ4cW2VWML95eIe7hOn9GZDfR+UtP428AeqQuGRbX4I87SgZ/3JTc+QsGArS2tnL//fcDmtPV2dlJW1sbhYV9K4Q/8YlPIKXk4YcfZsOGDQOWgzdt2pT6+sPEZDJxyimn8NFHH3Ho0CFmz57NmWeeSXV1NT/96U/p7u7mxRdfTHn+t771LZ588klAE7q12Wx4PB7y8vIyZuMgzAA+klJ6AIQQeUKIxVLKrUnG+oFdQoi36SsI/c2xMFQxfEKhEN3d3UZxRn/0qJ6UclwkV3RHNBAI4PF4iEQitLW1GRqBqc7Rq239fr+xr6qqGjAuEong9/uNbhN+v5/q6uOrLaQis2Sv86c3u/c3QMHJ2hIhpC+Me+ZDQ49RjC/CDKWXQf4I+mGe87uP3yOgRf6EKbtapdmnarmI/dj74S4WzJ2BAy/gYOfOnSxduhSn0zmg5/CePXtobm6muLg4qROVycgfwHnnnce7777Lb37zG9atWwdoS1aPP/54H3X7gwcP8uKLL/K9730PgHXr1lFZWckDDzzAvfdqP9JCodCoVC+m4L/Q+vbq+JIc01nTu2U9qtpXQ89hy83NxWKxGD1++/f6tdlsxGIxbLbRq6KPRCIcOHCAUChENBolEolQUVFBY2Mj69evx+l0EgqFaG1t5eDBg3R1dbFp0yYOHTpEa2srr732GvX19QNy83Sam5s5evQojz76KPF43Ihm9qeyshJA9cNVjIgJ4PzVa86fv15L6D/3cZh94/japsgMxWfDxW+M7Nyya/rdv1YrBMkm6ZQUOn+7d+/kSIObYDhGLObj/vvv55e//CVFRUXEYjGCwSAOh4Pm5mZuvvlm/vKXv/DNb36TdevWceWVfdvgZTLyB5rzd+utt3LXXXcxa9aslOMWLlxoOH4AJSUl3HLLLdx9990AtLe3U1JSMpZLTUImfJNKKeNCJBdRlFI+JYSwAQt7Dx2QUmblMoGq9u2L1Wod9/Zpfr+fjo4Oo/3XU089xd69e2lqamLr1q3k5eVhMpnw+Xx0d3cTCARobGykpqYGu91OaWkpZ555Ji6XC5fLZeTn6befeeYZGhoaOPHEE40xybY1a9ZQWVnJ2WefPa6vh2JiMjGcP9AqfaNesKWZoNu0DrZ+FS5+q2/eoCJ7OJYWZ95a6Nz1sdj3lLO0LZuQMa0IxeToo/W3e/durr96JeeuuIxIJMr3vvc9LrjgAgAuv/xyNm/ezPLly7n++ut56KGHOPnkk/nBD37At7/97QHOX7rcdNNNbNiwgba2NsrKyvjRj37E7bffPmDcokWLC7m74QAAIABJREFUsNvtfPvb3x7W/B9++CFLly417v/9739n1apVI7J1hFQLIb6JFu0DuBMtt28AQoiLgKeAWkAA5UKIL2ej1Mvxiu7HJ+b2ZRtz5szB5XJhsVg46aSTEEJwyy23MHv2bJxOJ93d3axbt46mpiYuu+wynn/+eVwuFytWrOD6669POa/T6WTGjBlJ/z8T0aOgCsVIyN5sV9csQHxc0Zmuxp+OxQWBRlX0kc28vRK2/MPIzj3yPGz6jFYwAnB0PQTcmbMtE0R9mhh1zNfn8O4P93DXbTewe9cu9u3bxy233GI8dtddd/HUU0/hcrl47733uOyyywAt9++9994bsSl//vOfaW5uJhKJ0NDQkPKL5de//jWrV68mJyfHONbe3s7Xv/51KioqWL16NQAHDhzg3//94zaKJSUl/P73v2f/fk1X+ZlnnuEf//EfR2zvCPg6sBxoBBqAc4E7Uox9CLhcSvlJKeUngCuAX42JlYo+RKNRo0K3f8cJvZdsOBwmGtVqd7K1l+y0adMoLi5m4cKFzJ49m2nTppGbmzvuUUqFIhXZ+7PBZIVz/vvjaI63BhDgOiG98/tHDhXZR0+lVsE9ElwJci85J8D6S+D0n8Pi4UWsRpUUWn9V1bUsWHBi0k41Z555JitXrjQqbseKqqoqVq1axQUXXMCXv/zlPo8VFxfz29/+ts+xk046ie9///vG/WuuuYZrrtGW4sPhMNdddx0LFy5krJBStgAp+v4NwCqlPJBw7kEhhCqFHAd0rb1EKRagTy6fvpdSZq3zp1BMNLLX+QOYn5Dm4q0ZXk6Xc6a2V1p/2UnUB0H38KtzdRKde13RY6RzjRaG89c3nayhZv/H7d2ScNttt42mVX249dZbKSwspLCw0EggP1ZsNtuYiEkfAzuEEL8H/tR7/2Zgxzjac1xjsViGXL5UTp9CkVmy2/nz1YOvDqZdCNM+MbwoUTI5EEX2oHdrGanDlpPg/MUCvXNlkcYfaJE9k2Wg3EuKQpDx4NZbbx1vE8aDbwB3Abq0yybgN+NnjkKhUIwt2e38VT4MVb+Dz3lg/leHf/7smyEvyxwChcZIdfl0nDMBoUV2I93HNtdo0l/rT0qtCMRsT7rsqxh9ejt7PNy7KRQKxXFHdjt/rjJteTDUrlVLplvpq3OW+mzPWhwzYN5twxd41jFZ4ZK/a9HgvT8B25S0I2pjmjvkKqNPXVU8At0faXqEjqljY0OWkUzbbKQIIe4Z4loPJ4x9Xkp5Y79WcIljT8uYYYpJhS7IHIlE6OzspL29nYaGBsxmM21tbUZBikIxUchy5693aa9lA2z+HCx/GuZ8cXhzxKPa0psiuyg5R9uOhemf1PaLvgXlqaUTEnE4HLS3t1NcXDw2DqA1v+99PQpoyiI9wjFESkl7e3smxZ+H0z7kn3v3Vw86KotQIs/jS1VVFXV1deTm5hIMBgmHw3R0dHDo0CFCoRAWi4X29nZsNht5eXljKWquUBwT2e0VGc5fr/xWupW+Ogf+E3b+M9zQ9XGLOEV2EOrQInXHsvTZ9r5WMTzvVshPr7K0rKyMhoYGWltbR37d4RCPacUdZqemaRj1QagNnDYtenkc4nA4KCsry8hcUsofDWNsc+/NO6WUfcrChRC/ALKoVFxDiTyPH/F4nA8//BC3282UKVOw2WyYzWZCoRBCCGw2Gzk5ObhcLmbMmIHD4VD9cxUThux2/nL6OX/pavzp2Kdo7bX0FnGK7OHN5VB4Olz47MjnqHsODv9WE1OedlFa+Z1Wq5W5c/9/e3ceH1dd73/89ZnJvjVJU5qudINCF8G2FBTZiqxSqYgLLleQC6I/lJ/3d1G4l3vx3iuK4gaiIiqLFxQUobIUoVigUGTRttDdltLSpk237Mkks31+f5yTOE2zTJKZOSeZz/PxmMfMnJw5552ZSeYz3/M93+8A30dD8c6D8JfPwIc2OYONr78F/n4TfLwNcgozl2OEE5EC4EpgNtDV/KKqPZ06fQ5HFnoX9LDMZLH29nbq6+spKChg+vTpnHDCCeTn5xMOh6moqGD+/PmUlZXx0ksvMW3atB67M8Tj8a7WwsbGRl8OVm2yk7+Lv4JxcNpjzoC+TZuhcNzAHl/kti5Y8ecvGnfO9p1w8dC2UzTJOXnitX+Gk3/pz5N7Es9KHnUctO6AgrFW+KXe/wKbcQZs/m+c4Vs2Ja4gIl/Emfljuoi8lfCjUmBVhnKaYaa8vJxp06Z1fWlsbGyksLCwa3iaWCzGvn37uub4jcVitLS0sHv3bl577TXq6uoIBoPs3bsXVaWsrKyv3RmTEf4u/gJBZ/quHQ9CyZSBHyK0gZ79qa3G6fs21LNzOwsr8N8wL526vwenfR7Gnu1dnpFrhqp+TEQudufu/Q3OEC6JfgM8DXwbuCFhebOq1mUqqBle4vE4bW1t7Nu3j9bWVmpra1m/fj11dXVEo9Gu+XzHjRtHWVkZFRUVFBcXM3bsWCZOnEhRURELFizgwIEDBINBZs+e7fWvZIzPiz+Ag687fcOOv37gj7WBnv2pxZ1ydajFX2FCvzE/DvMCUOhOU9hZ/I15n3MxqdY5knaDiMwBaoGjEldQ1UagUUS6H94tEZESVX03AznNMKGqbNu2jaamJkaNGsWOHTsACIVCNDY2MnHiRKqrqykvL6e6upp///d/Z/z48YgIoVCILVu2sHr1akKhENu3b6e2tpZJkybZgNXGF/xf/G3+AdT9DU7+xcAfG8yH478Go09KfS4zeJ3F31AP0ya2/BVOGNq20iWY5xzmbX3XOdy9dzlUnAiFY71ONtLcLSIVwE3A40AJ8J+9rPsUzlAvgtM/cCqwBae/oDHEYjEOHTpEY2MjpaWlzJ07l/nz51NSUkI4HGbdunW85z3voaioiNzcXJqamli9ejUrVqygsbGRcDhMU1MTGzZsoKGhgaqqKoLBIFOnTmXcuAF2XzImDfxf/BWMgZZtEG2DnKKBP/6930l9JjM0le+FOTf/45DoYBWMg/K5EG50ugj41elLobAaQrXwwvlw0k/hmC96nWpEUdVfujdXAn02A6vq3MT7IjIPpy+gySLxeJyWlhba29tpaWmhtbW16zoUCtHc3ExTUxPjx49n6tSpjB3rfGHrPqZfTU1N1xBS5eXlTJw4kVGjRhEOhykuLubAgQOcf/75PPXUU8ybN8+GgzG+4P/iL9rqXL9z/+A+MOMRd2gN+7blGxUnOpehCgSdgZ7bMzRsy2BVnexcH3jFuS6e4lmUkUpEvgV8V1Ub3PsVwP9T1Zv6e6yqrhaRk9Od0fiDqvLyyy+zadMmSkpKKCkpAZw5qQsLCyktLaWyspJQKERpaSnHHnss+fn5fW6vqqrqiKkS9+7dy+rVq6mrs+6kxn/8X/yNORW23wvlgxx8f/W/wDv/Cx9rSG0uM3gNG5xiPL9y6NvKH+1c/Kz+LTjw0j9mqLHiLx0uUNV/67yjqvUiciHOYeDDdJsVJADMA/akP6Lxg86BxouLi5kyZQqVlZUEg0FUlUgk0rVeOBwmEAhQWFhIVVWVh4mNST3/F3/TPg9jTkt6EN8jFE1y5n6NNNtAz37x5zNh4kfg5Lu9TpIZtc/CmuthlnueQfEABys3yQiKSL47by8iUgj01lyT+I8gitMH8A9pzmd8pPNwr4igqgQCAfLz8w+7hEIhRo0axVFHHWXDs5gRx//Fn8jgCz+wsf78JtLkHIb345h86dI1U81LkD8Gcoq9zTMyPQj8WUTude9fAdzf04qds4KISJlzV5szE9EhIkuADwFlwK9U9dlM7t/QNdjytGnTqK6u7nGd9vb2Hh8XCoXSms2YTPB/8TdUieOsWfHnvZZ3nGu/Ds2SDp3vwaM/AaPmeJtlhFLV77gDN3cOovg/qvpMT+uKyALgXtwWQBFpBK5U1b/2tx8RuQdnbuD9qjonYfn5wO1AEPilqt7aR9alwFK3X+L3ACv+PNLXdGyRSISOjg4OHjzI5s2baWlpIRQKdc3kYVO5meEsC4q/hJY/472Wt53rbCz+ggVQvcjbLCOYqj6NM4hzf+7Bmd/3JQAR+YC7LJmOxfcBdwK/7lwgIkHgJzjTxu0G3hCRx3EKwW93e/znVXW/e/sm93HGB5qamqirq6OlpYWWlhYOHTpEQ0MD+/bto7GxkZKSEqqqqrpOEikoKCAcDnsd25hBGfnFX+EEOOHbULnA6yQGEgZ4zqLDvoXjAIF133DmIB5KNwZzGBF5WVU/ICLNOGP3df0I55BuT521Yp2FH85KL4tItIf1jqCqK0VkSrfFC4FtqrrdzfQQcLGqfhunlbB7ZgFuBZ5W1dV9/G5XA1cDTJ5s/UTTbePGjUQiEYqKiqioqKC0tJSKigpmzpzJySfbyeBmZBn5xV8wD2bf0P96JjMmLHb6veWVe50kcwI58MEX4LkzoPY5K/5SSFU/4F73ezaXO54fwIsi8nPgtzgF4yeAF4YQYwKQOIfkbqCvauHLwAeBUSIyQ1Xv6mklVb0buBtgwYIF2tM6JnXi8Tjjx49n+nTni2lDQwN5eXldc/gaM5Jkx7s6tBfC9TBqltdJTNlM55JtxO0fZMO8pJx72HWDqh7Xz6rf73b/5oTbGSuuVPUO4I5M7c84w7uEw2FCoRBtbW0cOnSI0tJSm2rNZK3sKP7e+BI0b4UPrfc6ialZBqOOy64+fwDb3MadkimexhiJVDUmIltEZHJf8/Oq6llpilADJE5XM9FdNmQishhYPGPGjFRsLmvE43Ha2tpob2/vunSeqNE5rMuoUaMoLbXhv0x2yo7ir2gi7Hve6xQmHoOXlsBx/wonfsvrNJnV+f4rPtrbHCNXBbBBRF4HWjsXquqHO2+LyGdU9YFugzyTsO4PBrnvN4BjRGQqTtH3SeBTg9xW90xPAE8sWLDgqlRsbziKRCKHDb7cm3A4zMGDB6mtraW5ubnrcG1+fj5lZWUUFBRQUFBAIBCgrKyMoqIiAoFAuuMbn+r8MgDOXM5vvfUWGzdu9DBRZmVJ8WcDPftCaLcz3V62tfqBMw3d/pU2xl/6/EcS63Q++YP+JyAivwXOBKpEZDdws6r+SkSuBZ7BOcP3HlXdMNh9dNtf1rX8RaNRotEokUiEaDR62NAqiYdpo9EoHR0dXRdV5cCBAwSDQcaMGUNhYSH5+flHFHidY/yZ7KWq7Ny5s+uLxTHHHENjYyMdHR0ATJgwweOE6ZcdxV+p+4+zYR2Meb+3WbJZ15m+WVj8lU7ProGtM0xVX0xinZ+7/QObVPWHg9zPZb0sXwYsG8w2+9nfiG/5i8VixGKxrg/izmIvGAySn59PTk4Oubm5XYdyO4u9ziIuJyeHoqIi8vLyiMfjiAiVlSmYOtIMW4mtep327NnDypUrefHFF1m5ciW7d+/uGsh78eLFnH766Wzfvp2vfe1rWTGjS3YUf2NOd673PW/Fn5ea3TH+rAgyKTLQoV7c/oGXAYMq/kxqxONxotEojY2NBAIBcnJyCAQC5OXlkZub23W/04EDB7oO/Xafii0YDHatFwgEevzgNyPDQF/bzpbiD3/4w6xatQqA0aNHc9pppzFlyhS+853vEI1Guf322wG4664eT7wfkbKj+CuoglMfhqpTvE6S3Vq2g+RA4USvk5gRYiBDvSRYJSJ3Ag9zeP/AXsfc88pIOuzb2cIXj8dRVWKxGKpKYWEhRUVFhxVxiVSVSCRCYWEhpaWlNvSKGbBVq1Zx/vnnc+ONNzJr1iwCgQBr1qzhtttu8zqaZ7Lnr+joj3udwBx7LYw/HwI9/5M3Zijccfw+gNMC+LKqrull1RPd6/9OWKaA76ZfGc6HfTsLvM6CD5yWmEAg0FXohcPhI1rvetPZImhGnsG26A3E3LlzmTPHptfslD1/SZEWePdhGH0ylNsbwBNF452LMSkmIv8JfAx41F10n4j8XlW/2cPqV3bOxpHw+CzsiJoesViMaDR6WMGXk5NDMBg87FBuLBbzKqIZonQcWrcxFzMri85zj8PrX4B3f+d1kOy15U6of9PrFGZk+jRwkqrerKo3A6cAn+1l3Ud6WPb7tCUbAhFZLCJ3NzY2eh0laZ2Hc3Nzc7uGV8nNzbVhVbKUFXX+lD1/jbllzvy+tX/2Okl2CtfD374Mtcu9TmJGpj1AQcL9fLoNtCwix4nIR3GmVbsk4XJ5t8f6hqo+oapXjxo1yusoA9LZ2mcf/MNHJg69Gv/InsO+AGMXwabbbLw/L7S841xn4zAvJhMacQZ5Xo7Tf+8c4HURuQNAVb8CzAQuAsqBxQmPbQaGXZ86Y5JhRZ3pSXYVf9WLYOO34cDLMP4Cr9NklxZ3mBcr/kx6POZeOr3QfQVV/SPwRxF5n6r+JVPBjPEDK+pMouwq/qreD4F8aHjLir9My+YBnk3aqer9A1j9IyKyAQgBfwLeA3xVVR9IS7ghGElDvRhj/CN7+vwB5BTBJXth1te9TpJ9WrZDfpXT99KYFBORi0RkjYjUiUiTiDSLSFMvq5+rqk04h4B3ADOA6zOVdSCGa58/Y4y/ZVfxB5BX4XWC7DT/djjvDa9TmJHrR8DngNGqWqaqpd1n90iQ615/CPi9qg6fU2nNsDSQfnc2Q4nJhOwr/kK18OLFsOdpr5Nkl2ABlEzxOoUZuXYB6zW5T84nRGQzMB/4s4iMAdrTms6MaMn2p0v1esYMVnb1+QPIq4Ta56D4aOv3lynxCKz5ujPLik2xZ9Lja8AyEXkR6OhcqKo/6L6iqt4gIt8FGt25fluBizMX1fhdLBajpsYZKWgwhVhTUxOvv/46W7duJT8/n9zc3P4fZEwGZV/xF8yDo06DfSu8TpI92nbBlh9C+Vwr/ky63AK04IzXl9fTCiKySFVXiMglCcsSV3n0yEd5y074yCxVZf/+/dTU1NDQ0EBRURED6W8ZDof561//yrJly9i6dSujR4/mjDPOYMGCBfzxj39MY3JjBib7ij9wxvtb+3UI7YPCsV6nGfnsTF+TfuNVtb95G88AVnD4GH+dFB8Wf8N5bt/hpqWlhcbGRqqqqigtLaW6upry8vKk5h1WVTZs2MCyZctYu3Yt+fn5LFq0iIsuuogJEyZYPz7jO9lb/AHsex6mfNLbLNmg2cb4M2m3TETOVdVne1vBnfYNVb0ic7GM34VCIXbt2sWuXbuorq5mxowZVFRUkOyUeh0dHWzfvp0f/vCHhEIhZs+ezZIlS5g1a1aakxszeNlZ/FW81ykAgz0eHTKp1rIdAnlQON7rJGbk+iLwryISBsKAAJp4xq+I/EtfG+ipf6AZuWKxGDt37uTQoUMEg0HGjh3LzJkzqahIbkSI+vp6nnvuOWpqasjNzWXRokUsXryYk08+mZyc7PxoNcNHdr5DA0E42+b4zZiOA1Ay1XnejUkDVU1mvsbOdWYCJwGPu/cXA6+nI5fxH1WlsbGRQ4cOUV1dTXV1NRMmTKC9vT2pkzva29t57bXXWLZsGW+//TaqytFHH81NN91ESUlJBn4DY4YuO4u/TrF2UIWcQq+TjGyn3AOxjv7XM2aQxPnU/jQwVVX/R0QmAeNUtauoU9X/ctddCcxT1Wb3/jeApzKfeuRRVeLxeFJFVCQSSWqb7e3OKDz9bTMSiRCJRPo9s/aNN95g7969TJkyhTlz5lBcXNzjeuvXrwc4os/fk08+yapVq2hubmbevHls2bKFyZMnW+FnhpXsG+evU8s78PtyePdhr5Nkh2C+1wnMyPZT4H3Ap9z7LcBPell3LM6h4U5hd5nviMhiEbk72f5nXorH43R0OF/y+irAVJXW1lba2trIzc3tc92mpibq6+vJy8ujsLD3L+ltbW3s3r0bVaWysrLX/a5du5Y//elPFBUVccEFF/RY+MViMf7whz/w1FNPMXnyZObNm3fYz5ubmxk9ejSlpaVs3ryZYDDImWee2Ws2Y/woe1v+io92phqrXQHTLvc6zcgVrofXroaZX3GG2DEmPU5W1XkisgZAVetFpLdOvb8GXheRx9z7S4D7MpBxwIbL2b7RaJRIJEIgECAvL6/XVrp4PE5LSwvRaJTCwsJeC7p4PE59fT0dHR0UFxdTVlbW6zbr6+upq6sjPz+f6urqHvvbRSIRVq5cyaZNm6iqquKkk05i4sSJR6zX3NzMXXfdRU1NDQsXLmTx4sVHtPzFYjHWrl3Lzp07OfXUUznttNM4/vjj+3uKjPGV7C3+JABjz3LG+1MFG1E9PZq3wa5HYOpnvE5iRraIiARxhmzBnbUj3tOKqnqLiDwNdH4buUJV12Qm5sgTDoeJxWIEg0Hy8no/iS4SidDS0gJASUlJr+tGIhHq6uqIx+OUl5dTVFTU43rxeJz9+/fT2tpKaWkpY8aM6bFAbGpqYvny5ezdu5fZs2czZswYAoEjD3rt3buXVatWUVxczJIlS1i4cGGP23rqqafYtm0bM2bM4Oabb+b73/9+r7+zMX6VvcUfOGf8vvs7aN4KZcd6nWZksjH+TGbcATwGHCUitwCXAjf1trKqrgZWZyjbiKSqhMNh4vE4ubm5fZ7h2t7eTltbG8FgkJKSkl7HzguFQjQ0NBAIBBg9enSfBeLevXuJRCJUVVX1OhDzrl27eO655wiHw5x++unMmTOHbdu2EY1GD1vvlVde4YknnmDMmDFcddVVTJ48+Yhtvf322/zoRz+itraWWbNmcfrpp/faX9AYv8vu4q/6bOd63wor/tLFij+TAar6oIj8DTgbZ5iXJaq6yeNYI1YsFus6YSM/P7/HljRwCsS2tjY6OjrIy8ujuLi418O3jY2NtLa2kpeXR2VlZa/bbG1tZf/+/YgI48eP7/HQsary5ptv8uqrr1JYWMiFF17IpEmTjlgvGo3y2GOPsWbNGo466ig+85nP9Fj4vfjii9x7773k5ORw9tlnU1hY2Gs+Y4aD7C7+SqbDvB/AUWd6nWTkatkOBWMhx74hm/RS1c3AZq9zjHQD6d/X3NxMLBZLWf++uro66uvr++zfFw6HWblyJZs3b2bixIksWrSIsrKyI9ZraGjggQceYM+ePZxyyilUVlYe0ZIXjUZ54IEHePbZZ5kwYQLXXXcdy5cvZ8+ePf09Tcb4WnYXfyJw3Fe9TjGySQ5UzOt/PWOM7w2mf19paWmvZ/QOpH/fvn37aGtr67N/X2NjI8uXL6e2tpa5c+fyvve9r8ecu3bt4rHHHiMcDvPRj36U+fPn88orrxy2TkNDA3fccQcbN27kpJNO4pprrrHDvGbEyO7iDyAagtrnoOIEKD6yud8M0cKfeZ3AGDNEqkpHRweq2m//vo6ODlpbW/vt39fe3k59fT2BQICqqqpeC8RoNMqePXuIRqN99u/bt28fTz/9NB0dHZxxxhnMnj27x0Oz69at49lnn2XatGl87nOf6/Gs39raWr75zW/S0NDApZdeyiWXXGKHec2IYu/mjgOw8sOwy3dzuhtjspxfxvmLx+OoKnl5ef1OXdbe3k4wGKSsrKzXwg+cvnuBQIAxY8b0OdZfW1sbkUiEcePG9Vr4AWzdupVwOMx5553H3Llzey3W3nzzTUpKSrj22mu7Cr/9+/cTjUYpKCgA4NVXX6Wuro4vfOELXHrppVb4mRHH3tHFk6FkhnPSh0mtlu3wzCmwf6XXSYwZllT1CVW9uq+iJ5OSLYKCwWBSs3wEg8F+t6mqAH0eZu6Um5vbZ4HYub2SkpKuQ7ihUIitW7dSVlbG+PHju9YJBoNMmTKl330aMxxZ8QdQvQj2vwjxaP/rmuQ1bYVDr+GcfGmMMf6iqmzevBkR4bjjjkuqYDVmJLDiD5zx/iJNUG/jvKZUy9vOdcl0b3MYY0wP3nnnHZqbmzn22GO7Dvkakw2s+ANnqBcJOFORmdRp2Q7BAiis9jqJMcYcpq6ujt27dzN+/Hiqqqq8jmNMRlnxB1A4Ft7/Wxh3rnO/fb+3eUaKlu1QPNUprI0xxidisRhbtmyhuLiYadNsAHqTfexTudPRH3euD74KSyfDmzdBrN3bTMNd8ZR/zKJijDE+oKocPHiQWCzG8ccfb2fymqxk4/x1VzIdJn8cNtwCOx+GhT93TggxAzf/B14nMMaYwzQ2NtLR0cGMGTN6HVTamJHOvvJ0VzAG3v9rWLQcUFhxNrzxpczsu/5N2H4/7HseIs2Z2Wc6uUM0GGOMHzQ2NtLU1ERRURHV1dYX2WSvpFr+ROR84HYgCPxSVW9NayrX0jU13PbMFvY0hBhfXsj1581kyXsnpHwbPa/zQbhwHWz4JusO5HDNrSsItr3DRUet54Ljypg7NgeiLRBthVlfZ+nfA6x68T7OyX+SeE4Jx04Yz7RxYyGnFI65hqWbwvzmz89zVGQds0fVccGUNqbk74PWnfChDRDM4+1Xb2d6/b0AxDRAS8GxjJp0Bpz0MxA5Mue5x7BkbjnklkI8Aq9fTd2eteS0vs2ByCj26kQKpn+aBWde6xRioT0s3azc9uzfj3w+2vfDnmXQupOd725iT+12Gjpy+UP7p7norItZMlNh3wuQP4bnd8T58aoGDja3ESyawHXnzWXJxHdhx4MQ2g2tu6B+NWHN5bK3v0Vt3rxBvXapeG0zob8cfslpTDaLRqNs3ryZYDDI6NGjvY5jjKf6Lf5EJAj8BDgH2A28ISKPq+rGdAZbuqaGGx9dRygSA6CmIcSNj64DSPqDM5lt9LfO0viXuPHFdYQiIS4ctY2vVf4I9uNcAvmQU8zzkQ9x49PC+wrCTCqtpTgQouTg68QaOwhqB8vbzuLGJ1v4bPkK/u1op7g7dHAUdSXTqBy7AGKtLH3rAD96/SwkdgpH5+/lxKItzC/+O8cH1lK10Cn88l+7jFsrmjlYWs70/N3M2LiL3TXnMnHx4xDIpbHmNbbU5bO9/VQqgk1Mzd/DU2vWsHtUDUuOz4WlEzk3ns+xY8bTNrqACbn7+dFzVwANZez5AAAPb0lEQVT/lyWTdsOrVwBQGKmkSEdTVdBOy6FGbnx0HeM+uIuTa64G4CzgLPcl+Mi273HjozBp0RbmH3wQiidRG6liZf35vN1ezfrQdDraBv7apeK1zYR+3z8+yWlMttuyZQvhcJiqqiobz89kvWRa/hYC21R1O4CIPARcDKS1+LvtmS1dH5idQpEYtz2zJekPzWS20d86iT9/vmkBCzfeTyheQEVZBStvOAeAm25dQSgSYkVkISuaF3ZtZ0J5Iau+dhr/9d0XCUViPFq/iBea5rM7Mpa2eKHz80sWuTnWUBOqBOCd8AReaF7gbKOukFVuzs8UjOW00j1Mz9/Nto5JPHToPLYdmse33P1duPVOahpCRzwPE57ZwpLZ87it7joqYzuYml9DYaCDV1pOYGd7JS8/s4Ul178fFm/lzJ9sY0dDrNsWYtzw6kSe/+IWvvDLZ4i1H6Ay2Igi7AqPJRSN8ZXX5rDqBmeonI/euuKIHAN97fqTivdHJnL4JafxDxE5HrgOqAL+rKo2AXaa7dmzh0OHDjFt2jTy8/O9jmOM55Ip/iYAuxLu7wZO7r6SiFwNXA0wefLkIQfb00MR09fywW6jv3USfx7SAkJRZyDQloZwctsI5FLT0AHAwWgFB6MVg87xHS7nO7WXH/Zzga7ir89t5JXz093n0FMvPCHkjMdXOoOdDVt63MaO+jiUHcuz+7eiHDk0wp6G9oTbQ3/t+pOJfaQih19ymtQQkXuAi4D9qjonYXnSXWNUdRNwjYgEgF8DVvylUWtrK9u3b6eysrJrLl9jsl3KTvhQ1btVdYGqLhgzZsyQtze+vHBAywe7jf7WsW2kfhupkIl9pCKHX3KalLkPOD9xQULXmAuAWcBlIjJLROaKyJPdLke5j/kw8BSwLLPxs0ssFmPjxo3k5OQwc+ZMr+MY4xvJFH81wKSE+xPdZWl1/XkzKcwNHrasMDfI9ecl/weczDb6W8e2kfptpEIm9pGKHH7JaVJDVVcCdd0Wd3WNUdUw8BBwsaquU9WLul32u9t5XFUvAD6d2d8gu+zatYtQKMRxxx1Hbm6u13GM8Y1kDvu+ARwjIlNxir5PAp9Kayo4rE/eYM+STGYb/a1j20j9NlIhE/tIRQ6/5DRplVTXmE4iciZwCZBPHy1/qe5Kk43C4TD5+fmUl5d7HcUYX+m3+FPVqIhcCzyD05/lHlXdkPZkOB+cQ/2QTGYb/a1j20j9NlIhE/tIRQ6/5DT+oKovAC8ksd7dwN0ACxYssEEzjTEpk9Q4f6q6DOubYowxPUlb1xgRWQwsnjFjRio2Z4wxgM3wYYwxQ9XVNUZE8nC6xjyeig2r6hOqevWoUaNSsTljjAGs+DPGmKSJyG+BvwAzRWS3iFypqlGgs2vMJuB3qeoaIyKLReTuxsbGVGzOGGOAJA/7GmOMAVW9rJflaekao6pPAE8sWLDgqlRv2xiTvazlzxhjjDEmi4hq6k8iE5EDwM4BPKQKOJjyIKk3XHLC8MlqOVNvuGTtnvNoVR36CPEj0CD+pw6WH947fsgA/sjhhwxgOfyWAZLL0ev/1LQUfwMlIn9V1QVe5+jPcMkJwyer5Uy94ZJ1uOTMJn54TfyQwS85/JDBcvgvQypy2GFfY4wxxpgsYsWfMcYYY0wW8Uvxd7fXAZI0XHLC8MlqOVNvuGQdLjmziR9eEz9kAH/k8EMGsByJ/JABhpjDF33+jDHGGGNMZvil5c8YY4wxxmSAFX/GGGOMMVnE0+JPRM4XkS0isk1EbvAyS39EZIeIrBORtSLyV6/zdBKRe0Rkv4isT1hWKSLLRWSre13hZcZOvWT9hojUuM/rWhG50MuMbqZJIvK8iGwUkQ0icp273FfPax85ffWcikiBiLwuIm+6Of/LXT5VRF5z//4fdufFNcYYk2ae9fkTkSDwd+AcYDfO5OiXqepGTwL1Q0R2AAtU1Q+DO3YRkdOBFuDXqjrHXfZdoE5Vb3WL6gpV/bqXOd1cPWX9BtCiqt/zMlsiERkHjFPV1SJSCvwNWAJcjo+e1z5yfhwfPaciIkCxqraISC7wMnAd8C/Ao6r6kIjcBbypqj/zMqsxxmQDL1v+FgLbVHW7qoaBh4CLPcwzLKnqSqCu2+KLgfvd2/fjFASe6yWr76jqXlVd7d5uBjYBE/DZ89pHTl9RR4t7N9e9KLAIeMRd7vnzaY4kIktE5Bduy+y57rJpIvIrEXmkv8enMMcR+xSRySKy1D2ikPYjR71kOFNEXhKRu0TkzHRn6CPHEa9TpojILBH5nYj8TEQuzeS+EzIEROQWEfmxiHzOiwxujuPd98IjIvLFDO63WETud98Dn07mMV4WfxOAXQn3d+PDD64ECjwrIn8Tkau9DtOPsaq6171dC4z1MkwSrhWRt9x/4r44RN1JRKYA7wVew8fPa7ec4LPnVESCIrIW2A8sB94GGlQ16q7i97//Yaenbhbu8qS726jqUlW9CrgG+IS7bLuqXpnhHD3tcy7wiKp+Hue970UGxTmaUYDzHu5TunL09DolIxV5gAuAH6vqF4F/SnbfKc5wMTARiJDE65CuHKq6SVWvwTn6cupgcgwyzyU4fwtXAR9Oageq6skFuBT4ZcL9zwJ3epUnibwT3OujgDeB073OlJBtCrA+4X5Dt5/Xe52xj6xjgSDOF5FbgHu8zpiQrQTnUOolfn5ee8jp5+e0HHge+ABOy3/n8kmJ7wu7pOS5Ph2Y1+3vLYhTeE8D8tz/ZbNwCqknu12OSnjc94F53bb/iAc5Hkm4Pdp9L60ArvAoQ8C9Hgs86NVz0dfrlO487uUnwG3AKi/ep8ANwBcG8r5M12uDU3w9DXwqg3+/NwInuuv8Jpnt5+CdGpx/+J0must8SVVr3Ov9IvIYzmHrld6m6tU+ERmnqnvdfmH7vQ7UG1Xd13lbRH6B80fkObdv2h9w/qE/6i723fPaU06/PqcAqtogIs8D7wPKRSRHndY/X//9D0equtJtEU7U1d0GQEQeAi5W1W8DF3XfhogIcCvwtLpdDLzI0YsrgJvd7T8C3JvpDKoad2/WA/lJrJ+WHIN9nVKY5/+I04//0V5+ntYMIrIbCLt3YwPNkKoc7nYeBx4XkaeA3wwmy0Dz4LR2TgTWkuQRXS8P+74BHCPOGX95wCeBxz3M0yv3eHpp523gXGB934/y1ONAZ7+HzwF/9DBLn9wiqtNH8MHz6v4j/RWwSVV/kPAjXz2vveX023MqImNEpNy9XYhzktcmnFabzj5Cnj+fWWKg3W2+DHwQuFRErgEQkdHinKDzXhG5MRM5etnnn4CvuMt3eJFBRC4RkZ8D/wvcOYgMKclBD6/TEAw0zxQRuRv4NU7rXyoM9H36KHCeiPyY1DbKDPS5OFNE7nDfE8tSmKO/PI8CHxWRnwFPJLMhz1r+VDUqItcCz+A0Zd6jqhu8ytOPscBjzmctOTjNqn/yNpJDRH4LnAlUud9+bsb5Bvg7EbkS2InT/8BzvWQ9U0ROxOk7swP4gmcB/+FUnG4I68Tppwbwb/jvee0t52U+e07HAfe7LQMB4Heq+qSIbAQeEpFvAmtwClnjI6p6B3BHt2WHcPqWZTLHEftU1fX848uDVxkeZRCtXWnIccTrlME8OwBP+8GrahuQdD/UNOZ4AXjBg/224rSEJ83Lw76o6jLSUx2nlNvEeoLXOXqiqpf18qOzMxokCb1k9d0Hvqq+DEgvP/bN89pHTl/9TanqW/TQId/9u1qY+URZzS/dbfyQww8Z/JSjkx/y+CGDn3J0Slkem+HDGGOyh1+62/ghhx8y+CmHn/L4IYOfcqQ8jxV/xhgzArndLP4CzBSR3SJypXtyTWd3m004h+DT2t3GDzn8kMFPOfyUxw8Z/JQjU3k8m+HDGGOMMcZknrX8GWOMMcZkESv+jDHGGGOyiBV/ZshEpFxEvuTeHi8ZnPPTGGOMMQNjff7MkLmjkD+pqnM8jmKMMcaYfng6zp8ZMW4FprsDDW8FjlfVOSJyObAEKAaOAb6HMx/hZ4EO4EJVrROR6ThzQ44B2oCrVHVz5n8NY4wxZuSzw74mFW4A3lbVE4Hru/1sDnAJcBJwC9Cmqu/FOYX9n9x17ga+rKrzgX8FfpqR1MYYY0wWsuLPpNvzqtqsqgeARv4x7+A6YIqIlADvB37vthz+HGc6MGOMGdFE5JciMqufde4TkSOmsXPn1P1UinJMERHP51U3mWOHfU26dSTcjifcj+O8/wJAg9tqaIwxWUNV/3kID58CfAr4TWrSmGxiLX8mFZqB0sE8UFWbgHdE5GMA4vDlPMrGGNOdiFwvIl9xb/9QRFa4txeJyIPu7XNF5C8islpEfu8e8UBEXhCRBe7tK0Xk7yLyuoj8QkTuTNjN6SLyiohsT2gFvBU4TUTWishXu2V6SEQ+lHD/PhG51G3he8nNsVpE3t/D73N54r5F5EkRObOv38MMP1b8mSFT1UPAKvewwW2D2MSngStF5E1gA3BxKvMZY0wavQSc5t5eAJSISK67bKWIVAE3AR9U1XnAX4F/SdyAiIwH/gM4BTgVOK7bPsYBHwAuwin6wOlr/ZKqnqiqP+y2/sPAx91t5wFnA08B+4Fz3ByfAO5I9pdM5vcww4cd9jUpoapH9D1R1fuA+xLuT+npZ6r6DnB+ehMaY0xa/A2YLyJlON1aVuMUgacBX8Ep6GbhfEEGZ8SDv3TbxkLgRVWtAxCR3wPHJvx8qarGgY0iMjaJTE8Dt4tIPs7/1pWqGhKRUcCdInIiEOu2j/4k83uYYcKKP2OMMWaQVDUiIu8AlwOvAG8BZwEzgE3AdGC5ql42hN0k9p2WJDK1i8gLwHk4LXwPuT/6KrAPOAHnyF97Dw+PcvhRwYKE/Q719zA+YYd9jTHGmKF5CWeYqpXu7WuANerMovAqcKqIzAAQkWIR6d7i9gZwhohUiEgO8NEk9tlfX+uHgStwWiD/5C4bBex1WxE/CwR7eNwO4EQRCYjIJJxWSZL8PcwwYcWfMcYYMzQv4fTL+4uq7sNpUXsJwB3m6nLgtyLyFs6h0sP69KlqDfAt4HVgFU4B1tjPPt8CYiLyZvcTPlzPAmcAz6lq2F32U+Bzbv/q44DWHh63CngH2IjTJ3B1sr+HGT5sejdjjDHGYyJSoqotbsvfY8A9qvqY17nMyGQtf8YYY4z3vuEOdL8ep+Vtqcd5zAhmLX/GGGOMMVnEWv6MMcYYY7KIFX/GGGOMMVnEij9jjDHGmCxixZ8xxhhjTBax4s8YY4wxJov8f11VdT6oYfXDAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 648x216 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# make fancy plots\n",
"times = np.arange(n_observations)\n",
"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9,3))\n",
"# plot data and time trace of one of the posterior probabilities\n",
"ax1.scatter(times, observations)\n",
"ax1.plot(times, target_samples.sum(1) / n_samples,\n",
" ls='--', color='orange', label=r'$p(x=1|y_{1:t})$')\n",
"ax1.set_yticks((0, 1))\n",
"ax1.set_xlabel(\"time\")\n",
"ax1.legend()\n",
"\n",
"# for each time step, plot cumulative distribution of weights\n",
"# whose lines get darker with increasing time step\n",
"for alpha, weights in zip(np.linspace(1, 0, n_observations),\n",
" all_weights):\n",
" ax2.plot(np.sort(weights),\n",
" np.linspace(0, 1, n_particles, endpoint=False),\n",
" color=\"black\", alpha=alpha)\n",
"# double log scale because the plot sucks otherwise\n",
"ax2.set_xscale(\"log\")\n",
"ax2.set_yscale(\"log\")\n",
"ax2.set_ylabel(\"empirical cumulative\\ndistribution function\")\n",
"ax2.set_xlabel(\"weight value\")\n",
"\n",
"fig.tight_layout()\n",
"plt.show()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The left plot shows the data together with the posterior probabilities for the hidden state to be $1$. The right plot shows the empirical cumulative distribution function (eCDF) of the weight values for different times. So the $y$-axis shows the fraction of weights which are smaller than a certain value. The more time progresses, the fainter the lines get. They shift towards the left, meaning that more and more weights have values of essentially zero, which nicely illustrates the weight decay."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "venv",
"language": "python",
"name": "venv"
},
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Raw
smc.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import logsumexp
choice = np.random.choice
class Pdf(object):
def __call__(self, x):
"""log-probability"""
pass
def sample(self, n):
pass
class ImportanceSampler(object):
def __init__(self, posterior):
self.posterior = posterior
def sample(self, imp_dist, n_samples):
imp_samples = imp_dist.sample(n_samples)
imp_weights = self.calculate_weights(imp_dist, imp_samples)
# doing stuff in log-space to avoid underflows
normalized_imp_weights = imp_weights - logsumexp(imp_weights)
return imp_samples, normalized_imp_weights
def calculate_weights(self, imp_dist, imp_samples):
return self.posterior(imp_samples) - imp_dist(imp_samples)
class SequentialImportanceSampler(ImportanceSampler):
def __init__(self, posterior):
super(SequentialImportanceSampler, self).__init__(posterior)
self.weights = None
def sample(self, n_samples):
(imp_samples, weights) = super(
SequentialImportanceSampler, self).sample(self.posterior.prior,
n_samples)
# store current weights to make them available to recursion
self.weights = weights
return imp_samples, weights
def calculate_weights(self, _, imp_samples):
if self.weights is None:
self.weights = np.log(np.ones(len(imp_samples)) / len(imp_samples))
likelihoods = self.posterior.likelihood_function(imp_samples)
unnormalized_weights = self.weights + likelihoods
self.weights = unnormalized_weights - logsumexp(unnormalized_weights)
return self.weights
class Gaussian(Pdf):
def __init__(self, mu=0, sigma=1):
self.mu = mu
self.sigma = sigma
def __call__(self, x):
return -0.5 * (x - self.mu) ** 2 / self.sigma / self.sigma
def sample(self, n):
return np.random.normal(self.mu, self.sigma, n)
class Uniform(Pdf):
def __init__(self, low, high):
self.low = low
self.high = high
def __call__(self, x):
return np.repeat(-np.log(self.high - self.low), len(x))
def sample(self, n):
return np.random.uniform(self.low, self.high, n)
class HMMPosterior(Pdf):
def __init__(self, trans_matrix, obs_probs):
self.observations = []
self.obs_probs = obs_probs
self.prior = trans_matrix
def likelihood_function(self, hidden_states):
return np.log(self.obs_probs[hidden_states, self.observations[-1]])
def add_observation(self, x):
self.observations.append(x)
class TransitionMatrix(object):
def __init__(self, transition_matrix):
self.transition_matrix = transition_matrix
self.hidden_states = np.arange(len(transition_matrix))
def sample(self, old_states):
return np.array([choice(self.hidden_states,
p=self.transition_matrix[hidden_state])
for hidden_state in old_states])
class SequentialPrior(Pdf):
def __init__(self, distribution, first_state):
self.distribution = distribution
self.samples = first_state
def sample(self, n_samples):
self.samples = self.distribution.sample(self.samples)
return self.samples
if True:
if not False:
# importance sampling
n_samples = 200000
target = Gaussian()
imp = Uniform(-10, 10)
sampler = ImportanceSampler(target)
biased_samples, logws = sampler.sample(imp, n_samples)
# we have samples from the importance distribution,
# but not from the target distribution. So we have
# to reweight. To get samples from the target,
# just draw from the list of importance samples with
# probabilities given by the weights
samples = choice(
biased_samples, 2*n_samples, p=np.exp(logws))
fig, ax = plt.subplots()
ax.hist(samples, bins=30)
plt.show()
if True:
# sequential importance sampling
n_particles = 100
# define multinomial HMM (discrete hidden states,
# discrete observation states)
hidden_states = [0, 1]
observation_states = [0, 1]
start_prob = np.array([0.5, 0.5])
trans_matrix = np.array([[0.9, 0.1],
[0.1, 0.9]])
obs_probs = np.array([[0.9, 0.1],
[0.1, 0.9]])
# draw first state of Markov chain from the start distribution
start_state = choice(hidden_states, n_particles, p=start_prob)
pdist = TransitionMatrix(trans_matrix)
prior = SequentialPrior(pdist, start_state)
posterior = HMMPosterior(prior, obs_probs)
sampler = SequentialImportanceSampler(posterior)
# generate data
n_observations = 30
# we have to start with some state...
true_states = [0]
observations = []
for _ in range(n_observations):
new_state = choice(hidden_states, p=trans_matrix[true_states[-1]])
new_observation = choice(
observation_states, p=obs_probs[new_state])
observations.append(new_observation)
true_states.append(new_state)
# will store sucessive posterior samples; is enlarged in each iteration
states_samples = None
# stores log-weights for all samples
all_log_ws = []
for observation in observations:
# we're on-line: a new observation is arriving!
posterior.add_observation(observation)
# draw samples from importance distribution and
# store them along with their corresponding importance
# weights
(samples, log_ws) = sampler.sample(n_particles)
if states_samples is None:
states_samples = samples.reshape(1, n_particles)
else:
# append n_particles samples to states_samples
states_samples = np.vstack((states_samples, samples))
all_log_ws.append(log_ws)
# ... and back to actual probabilities
all_weights = np.exp(all_log_ws)
# get samples from target by reweighting importance samples
target_samples = np.array([choice(samples, n_samples, p=weight)
for weight, samples in zip(all_weights,
states_samples)])
# make fancy plots
times = np.arange(n_observations)
fig, (ax1, ax2) = plt.subplots(1, 2)
# plot data and time trace of one of the posterior probabilities
ax1.scatter(times, observations)
ax1.plot(times, target_samples.sum(1) / n_samples,
ls='--', color='orange', label=r'$p(x=1|y_{1:t})$')
ax1.set_yticks((0, 1))
ax1.set_xlabel("time")
ax1.legend()
# for each time step, plot cumulative distribution of weights
# whose lines get darker with increasing time step
for alpha, weights in zip(np.linspace(1, 0, n_observations),
all_weights):
ax2.plot(np.sort(weights),
np.linspace(0, 1, n_particles, endpoint=False),
color="black", alpha=alpha)
# double log scale because the plot sucks otherwise
ax2.set_xscale("log")
ax2.set_yscale("log")
ax2.set_ylabel("empirical cumulative\ndistribution function")
ax2.set_xlabel("weight value")
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment