{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Particle Filter Resamplers: Tutorial\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\nStone Soup comes with a number of resamplers, for the Particle Filter, that can be\nused straight out of the box. This example explains each of the resamplers and compares the\nresults of using each one.\n\n\n**Resamplers currently available in Stone Soup:**\n\n- :class:`~.SystematicResampler`\n- :class:`~.MultinomialResampler`\n- :class:`~.StratifiedResampler`\n- :class:`~.ResidualResampler`\n- :class:`~.ESSResampler` (Effective Sample Size Resampler)\n\nThe last two resamplers (Residual and ESS) are preprocessing methods that require the use of\nanother resampler.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plotter for this notebook\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n\n\ndef plot(normalised_weights, u_j=None, stratified=False, residual=False):\n nparts = len(normalised_weights)\n floors = np.floor(normalised_weights * nparts)\n\n # Plot\n fig, ax = plt.subplots()\n fig.set_size_inches(10, 1.5)\n plt.xlim([0, 1])\n for i, particle_weight in enumerate(normalised_weights):\n l = 0\n for j in range(i):\n l += normalised_weights[j]\n ax.barh(['CDF'], [particle_weight], left=l)\n if residual is not False:\n if floors[i] != 0:\n ax.barh(['CDF'], [(1 / nparts) * floors[i]], left=l, color='White')\n\n if u_j is not None:\n plt.plot([u_j, u_j], [-0.4, 0.4], 'k-', lw=4)\n if stratified is not False:\n plt.plot([s_lb, s_lb], [-0.4, 0.4], 'w:', lw=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cumulative Distribution Function\n\nThe Systematic, Multinomial, and Stratified resamplers all use a cumulative distribution\nfunction (CDF) to randomly pick points from. The CDF is calculated from the cumulative sum\nof the normalised weights of the existing particles.\n\n**Example**\n\nGiven the following set of particles and their weights:\n\n.. list-table:: Particle Weights\n :widths: 25 25 50\n :header-rows: 1\n\n * - Particle\n - Weight\n - Normalised Weight\n * - 1\n - 1\n - 0.1\n * - 2\n - 5\n - 0.5\n * - 3\n - 1.5\n - 0.15\n * - 4\n - 0.5\n - 0.05\n * - 5\n - 2\n - 0.2\n\nWe would calculate the following CDF:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"particle_weights = [1, 5, 1.5, 0.5, 2]\nn_particles = len(particle_weights)\n\nnormalised_weights = np.array([particle/sum(particle_weights) for particle in particle_weights])\nplot(normalised_weights)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Each of the resamplers then use a different method to pick points in the CDF. For example,\nassume the following points are picked: 0.15, 0.41, 0.57, 0.63, and 0.89. These points are\nshown on the CDF below.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"u_j = [0.15, 0.41, 0.57, 0.63, 0.89]\nplot(normalised_weights, u_j)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As shown above, the points 0.15, 0.41, and 0.57 all fall within the section of the CDF\ncorresponding to particle 2, while the points 0.63 and 0.89 fall within the sections of\nthe CDF corresponding to particles 3 and 5, respectively.\n\nHence, for this example, the number of times each particle is resampled is as follows:\n\n.. list-table:: Particle Weights\n :widths: 25 25 50\n :header-rows: 1\n\n * - Particle\n - No. of times resampled\n - Weight of new resampled particles\n * - 1\n - 0\n - 0.2\n * - 2\n - 3\n - 0.2\n * - 3\n - 1\n - 0.2\n * - 4\n - 0\n - 0.2\n * - 5\n - 1\n - 0.2\n\nAs shown in the table, the resampler assigns a new weight to each sample - by default, each of\nthe resamplers included in Stone Soup give all particles an equal weight of $1/N$ where\n$N$ = no. of resampled particles.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Resamplers\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Multinomial Resampler\n\nThe :class:`~.MultinomialResampler` calculates $N$ independent random numbers from the\nuniform distribution $u \\sim U(0,1]$, where $N$ is the target number of particles to\nbe resampled. In most cases, we use $N = M$, where $M$ is the number of existing\nparticles to be resampled. However, the Multinomial resampler can upsample ($N > M$) or\ndownsample ($N < M$) to a value $N \\neq M \\in \\mathbb{N}$. The Multinomial resampler\nhas a computational complexity of $O(MN)$ where $N$ and $M$ are the number of\nresampled and existing particles respectively.\n\nAn example of how the Multinomial resampler picks points from the CDF is shown below. Black\nlines represent the chosen points.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n\n# Pick N random points in the uniform distribution u~U(0,1]\nu_j = np.random.rand(5)\n\nplot(normalised_weights, u_j)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Systematic Resampler\n\nUnlike the Multinomial resampler, the :class:`~.SystematicResampler` doesn't calculate all points\nindependently. Instead, a single random starting point in the range $[0,1/N]$ is chosen. $N$ points are then\ncalculated at equidistant intervals along the CDF, so that there is a gap of $1/N$\nbetween any two consecutive points. The Systematic resampler has a computational complexity\nof $O(N)$ where $N$ is the number of resampled particles.\n\nAn example of how the Systematic resampler picks points from the CDF is shown below. Black\nlines represent the chosen points.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Pick a starting point\ns = np.random.uniform(0, 1/n_particles)\n\n# Calculate N equidistant points from the starting point\nu_j = s + (1 / n_particles) * np.arange(n_particles)\n\nplot(normalised_weights, u_j)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Stratified Resampler\n\nThe :class:`~.StratifiedResampler` splits the whole CDF into $N$ evenly sized strata\n(subpopulations). A random point is then chosen, independently, from each stratum. This results\nin the gap between any two consecutive points being in the range $(0, 2/N]$. The\nStratified resampler has a computational complexity of $O(N)$ where $N$ is the\nnumber of resampled particles.\n\nAn example of how the Stratified resampler picks points from the CDF is shown below. Black\nlines represent the chosen points, with white dashed lines representing the strata boundaries.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Calculate lower bound of each stratum\ns_lb = np.arange(n_particles) * (1 / n_particles)\n\n# Independently pick a point from each stratum\nu_j = np.random.uniform(s_lb, s_lb + (1 / n_particles))\n\nplot(normalised_weights, u_j, s_lb)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Preprocessing methods\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Residual Resampler\n\nThe :class:`~.ResidualResampler` consists of two stages.\n\n**Stage 1**\n\nThe first stage determines which particles have weight $w^{i} \\geq 1/N$, where\n$i \\in 1, ..., N$ denotes each particle. Each of these particles is then resampled\n$N^{i}_{j} = floor(Nw^{i}_{j})$ times, where $j \\in 1, 2$ denotes the stage. Hence,\n$N_1 = \\sum_{i=1}^{N}N^i_1$ represents the number of particles that are sampled in stage 1.\nAs these weights have been represented in the resampled set of particles, we are only interested\nin the residual weights, left after the floor weights ($N^{i}_{1}$) have been subtracted.\n\n\nReconsider our example with 5 particles, with weights shown by the plot below.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plot(normalised_weights)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we're interested in particles of weight $w^{i} \\geq 1/N$, in this example we're looking\nfor particles with $w^{i} \\geq 1/5$. Hence, Particle 2 (orange, weight = 0.5) and Particle\n5 (purple, weight = 0.2). We resample these particles $N^{i}_{j} = floor(Nw^{i}_{j})$\ntimes:\n\n- Particle 2 gets resampled $N^{2}_{1} = floor(Nw^{2}_{1}) = floor(5 * 0.5) = 2$ times.\n- Particle 5 gets resampled $N^{5}_{1} = floor(Nw^{5}_{1}) = floor(5 * 0.2) = 1$ times.\n- Particles 1, 3 and 4 do not get resampled as their weights are $\\leq 1/5$.\n\n\nA total of 3 particles were sampled from stage 1 (2x Particle 2, 1x Particle 5), hence\n$N_1 = \\sum_{i=1}^{N}N^i_1 = 3$\n\n\nOriginal CDF with floor weights removed:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plot(normalised_weights, residual=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"New CDF:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"floors = np.floor(normalised_weights * n_particles)\nresiduals = (normalised_weights - floors/n_particles)\nplot(residuals)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The CDF of the residual weights is carried on to stage 2.\n\n\n**Stage 2**\n\n\nWe now look to resample the remaining $N_2 = N - N_1$ particles during stage 2. (Reminder:\nN is the total number of particles we want to sample, while $N_1$ is the number of\nparticles we have already resampled from stage 1). The residual weights from each particle are\ncarried over from stage 1. These residual weights are normalised, and used to calculate a CDF,\nwhich is shown below.\n\n\nNormalised residual weight CDF:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"normalised_residual_weights = residuals*(1/sum(residuals))\nplot(normalised_residual_weights)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We then use either the :class:`~.MultinomialResampler`, :class:`~.StratifiedResampler`, or\n:class:`~.SystematicResampler` to sample the remaining $N_2$ particles using the CDF.\n\n\nContinuing with our example, we want to resample the remaining $N_2 = N - N_1 = 5 - 3 = 2$\nparticles in stage 2. Here we will choose to use the :class:`~.MultinomialResampler` (we could\nalso use :class:`~.StratifiedResampler` or :class:`~.SystematicResampler`) to sample\nthe $N_2 = 2$ particles from the normalised residual weight CDF - more detail on how this\nis done can be found in the Multinomial Resampler section above.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"u_j = np.random.rand(2)\n\nplot(normalised_residual_weights, u_j)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This method reduces the number of particles that are sampled through the more computationally\nexpensive methods seen above. The Residual resampler also guarantees that any particle of weight\ngreater than $1/N$ will be represented in the set of resampled particles.\n\nWhen using the Residual resampler in Stone Soup, the Resampler requires a property\n'residual_resampler'. This property defines which resampler method to use for resampling the\nresiduals in stage 2.\nThe variable must be a string value from the following: 'multinomial', 'systematic',\n'stratified'. If no `residual_method` variable is provided, the multinomial method will be used\nby default.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### ESS Resampler\n\nThe :class:`~.ESSResampler` (Effective Sample Size) is a wrapper around another resampler. It\nperforms a check at each time step to determine whether it is necessary to resample the\nparticles. Resampling is only performed at a given time step if a defined criterion is met.\nKish's ESS criterion is used, which states that:\n\n\\begin{align}ESS = \\left(\\sum_{i=1}^{N} (W_{n}^i)^2\\right)^{-1}\\end{align}\n\nis less than or equal to the threshold, a float which is set by default to $\\frac{N}{2}$.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example in Stone Soup\n\n### Generate some particles\n\nAn example of resampling particles using the Stone Soup resamplers. We generate some particles\nusing the :class:`~.Particle` class. In this example, we give every particle an equal\nweight - each particle will have the same likelihood of being resampled. The resample function\nreturns a :class:`~.ParticleState`. The state_vector property of which contains a list of\nparticles to be resampled.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from stonesoup.types.particle import Particle\n\nparticles = [Particle(np.array([[i]]), weight=1/5) for i in range(5)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Simple Example\n\nIn this example, we use the :class:`~.MultinomialResampler` to resample the particles generated\nabove. We simply define the Resampler and call its `resample` method.\n\nThere are situations where we may want to resample to a different number of particles - for\nexample if you want more granular information, or want to be more computationally efficient.\nExcluding the :class:`~.ResidualResampler`, all Stone Soup resamplers can up- or down-sample as\nshown below.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from stonesoup.resampler.particle import MultinomialResampler\n\n# Resampler\nresampler = MultinomialResampler()\n\n# Resample particles\nresampled_particles = resampler.resample(particles)\nprint(\"---- State vector of resampled particles ----\")\nprint(resampled_particles.state_vector)\nprint(\"------ Weights of resampled particles -------\")\nprint(resampled_particles.weight)\n\n# Repeat while upsampling particles\nupsampled_particles = resampler.resample(particles, nparts=10)\nprint(\"---- State vector of upsampled particles ----\")\nprint(upsampled_particles.state_vector)\nprint(\"------ Weights of upsampled particles -------\")\nprint(upsampled_particles.weight)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we can see which particles were chosen to be resampled and the weights of the new\nparticles.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example using ESS and ResidualResampler\n\nIn this example, we use both the Effective Sample Size method, and the residual resampler, using\nthe systematic method to resample the residuals.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from stonesoup.resampler.particle import ESSResampler, ResidualResampler\n\n# Define Resampler\nsubresampler = ResidualResampler(residual_method='systematic')\nresampler = ESSResampler(resampler=subresampler)\n\n# Resample Particles\nresampler.resample(particles)"
]
}
],
"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.12.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}