Particle filters: How to do resampling?

Question

I understand the basic principle of a particle filter and tried to implement one. However, I got hung up on the resampling part.

Theoretically speaking, it is quite simple: From the old (and weighted) set of particles, draw a new set of particles with replacement. While doing so, favor those particles that have high weights. Particles with high weights get drawn more often and particles with low weights less often. Perhaps only once or not at all. After resampling, all weights get assigned the same weight.

My first idea on how to implement this was essentially this:

1. Normalize the weights
2. Multiply each weight by the total number of particles
3. Round those scaled weights to the nearest integer (e.g. with int() in Python)

Now I should know how often to draw each particle, but due to the roundoff errors, I end up having less particles than before the resampling step.

The Question: How do I "fill up" the missing particles in order to get to the same number of particles as before the resampling step? Or, in case I am completely off track here, how do I resample correctly?

Answer 1

The issue you're running into is often referred to as sample impoverishment. We can see why your approach suffers from it with a fairly simple example. Let's say you have 3 particles and their normalized weights are 0.1, 0.1, 0.8. Then multiplying by each weight by the 3 yields 0.3, 0.3, and 2.4. Then rounding yields 0, 0, 2. This means you would not pick the first two particles and the last one would be picked twice. Now you are down to two particles. I suspect this is what you have been seeing when you say "due to the roundoff errors, I end up having less particles."

An alternative selection method would be as follows.

1. Normalize weights.
2. Calculate an array of the cumulative sum of the weights.
3. Randomly generate a number & determine which range in that cumulative weight array to which the number belongs.
4. The index of that range would correspond to the particle that should be created.
5. Repeat until you have the desired number of samples.

So, using the example above we would start with the normalized weights. We would then calculate the array [0.1, 0.2, 1]. From there we calculate 3 random numbers say 0.15, 0.38, and 0.54. This would have us pick the second particle once and the third particle twice. The point is that it gives the smaller particles a chance to propagate.

One thing to note is that while this method will deal with impoverishment it can lead to a suboptimal solutions. For instance, it may be that none of the particles really match your given location well (assuming you're using this for localization). The weights only tell you which particles match best, not the quality of the match. As such when you take additional readings and repeat the process you may find that all your particles group at a single location that is not the correct location. This is usually because there were no good particles to start.

Answer 2

As I guess you found out yourself, the resampling method you are proposing is slightly flawed, as it should not alter the number of particles (unless you want to). The principle is that the weight represents the relative probability with respect to the other particles. In the resampling step, you draw from the set of particles such that for each particle, the normalized weight times the number of particles represents the number of times that particle is drawn on average. In that your idea is correct. Only by using rounding instead of sampling, you will always eliminate particles for which the expected value is less than half.

There are a number of ways to perform the resampling properly. There is a nice paper called On resampling algorithms for particle filters, comparing the different methods. Just to give a quick overview:

• Multinomial resampling: imagine a strip of paper where each particle has a section, where the length is proportional to its weight. Randomly pick a location on the strip N times, and pick the particle associated with the section.

• Residual resampling: this approach tries to reduce the variance of the sampling, by first allocating each particle their integer floor of the expected value, and leave the rest to multinomial resampling. E.g. a particle with an expected value of 2.5 will have 2 copies in the resampled set and another one with an expected value of 0.5.

• Systematic resampling: take a ruler with regular spaced marks, such that N marks are the same length as your strip of paper. Randomly place the ruler next to your strip. Take the particles at the marks.

• Stratified resampling: same as systematic resampling, except that the marks on the ruler are not evenly placed, but added as N random processes sampling from the interval 0..1/N.

So to answer your question: what you have implemented could be extended to a form of residual sampling. You fill up the missing slots by sampling based on a multinonmial distribution of the reminders.

Answer 3

For an example of python code that properly implements resampling, you might find this github project to be useful: https://github.com/mjl/particle_filter_demo

Plus, it comes with its own visual representation of the resampling process, that should help you debug your own implementation.

In this visualization, the green turtle shows the actual position, the large grey dot shows the estimated position and turns green when it converges. The weight goes from likely (red) to unlikely (blue).

Answer 4

one simple way to do this is numpy.random.choice(N, N, p=w, replace=True) where N is the no. of particles and w = normalized weights.

Answer 5

I use @narayan's approach to implement my particle filter:

new_sample = numpy.random.choice(a=particles, size=number_of_particles, replace=True, p=importance_weights)

a is the vector of your particles to sample, size is the count of particles and p is the vector of their normalized weights. replace=True handles bootstrap sampling with replacement. The return value is a vector of new particle objects.