Skip to content

Particle Filter

Jump to bottom
Bob Stienen edited this page Jun 5, 2020 · 3 revisions

Particle filter implementation documentation

The implementation of the particle filter algorithm in the high dimensional sampling (HDS) framework is highly configurable, in order to get optimal performance a dedicated search for the best parameters might therefore be needed. In this document the workings of the algorithm are explained. Configurable parts of the algorithm are discussed where they influence its workings.

Defining the test function

As the algorithm is implemented within the HDS framework, it expects the functions to sample to be child classes of the functions.TestFunction class. Writing such a class is relatively easy. The following code implements a 3-dimensional quadratic function for example:

class ExampleFunction(func.TestFunction):
    def __init__(self):
        self.ranges = [[-10, 10], [-10, 10], [-10, 10]]
        super(ExampleFunction, self).__init__()

    def _evaluate(self, x):
        return np.sum(x**2, axis=1)

    def _derivative(self, x):
        raise func.NoDerivativeError()

You can write such a function (i.e. class) yourself by following the steps outlined in the manual.

Initial sampling (sampling the seed)

In order for the particle filter to work, it needs an initial set of samples. The HDS implementation automatically performs this initial sampling the first time the algorithm is run. It needs two pieces of information for this, all of which are defined at construction of the algorithm class:

  • seed_size: the number of samples to take;
  • seed_distributions: a list of strings, 'log' or 'linear', indicating for each dimension in the problem which distribution should be used for sampling;

The test function which needs to be sampled is not provided at construction, but instead when the algorithm is called.

Sampling ranges

The ranges within which the algorithm will sample are defined through the ranges parameter at construction of the algorithm class. For example:

ranges = [[-3, 3], [51, 102]]

defines a 2-dimensional plane bounded by -3 and 3 for the first dimension, and 51 and 102 for the second dimension. If ranges is set to None (which is the default), the ranges defined in the test function will be used. If one of the dimensions is unbounded (indicated by np.inf or -np.inf in the ranges), the maximum (minimum) is set to the value of inf_replace (-inf_replace), which can be provided during construction.

During the initial sampling these boundaries are always respected. Subsequent sampling of values might go outside of these boundaries however. If this is not allowed, the parameter hard_ranges can be set to True during construction. This will guarantee that all samples will lay within the boundaries, during any iteration.

Custom samplings

Instead of having the particle filter perform the initial sampling itself, the seed can also be defined by the user themselves:

ParticleFilter.set_seed(data, function)

As can be seed, the function which will be sampled needs to be known.

Particle filter sampling

When an initial data set is known to the algorithm, all subsequent calls to the class use the particle filter algorithm. The number of samples taken in each of these calls is defined by iteration_size.

During a call, new samples are taken from gaussians drawn over the samples from the previous iteration. The standard deviation of these gaussians can be configures by the user in the form of a gaussian constructor.

Gaussian constructor

A gaussian constructor is a function of the form

standard_deviations = gaussian_constructor(algorithm, previous_samples, function_values)

There are two of these functions readily available within the particlefilter module: gaussian_constructor_linear (default) and gaussian_constructor_log. The first one of these multiplies a width parameter w with the total range for each of the dimensions and takes the result as standard deviations, whereas the latter multiplies the result by the values for the samples.

Width parameter w

The width of the gaussians can be configured such that it automatically decays over time. At construction the initial value is defined by width, the decay of this width is given by width_decay. How these two interact is defined by the width_scheduler. This is a function that takes two arguments and returns the new width:

new_width = scheduler(algorithm, width_decay)

In this function, any update of the width should also be written to algorithm.width.

There are two prebaked functions like this available in the particlefilter module: width_schedule_exponential (default) and width_schedule_exponential_10stepped. The first one performs the following operation after each iteration of the algorithm:

algorithm.width = algorith.width * width_decay

The 10-stepped version does the same, but only after each 10 iterations.

Sample selection

In order for the algorithm to approach the minimum of the test function, all samples (i.e. 'particles') from the previous iteration are subjected to a selection function. This selection influences which samples are used as means for the sampling in the current iteration. In other words: not all samples from the previous iteration are used, a selection is made.

Selection methods can be defined by the user and provided as selector_function during construction of the algorithm class. Such a function should have the following form:

selected_indices, previous_samples, previous_function_values = selector_function(algorithm, select_n_samples, samples, function_values)

Four of these functions are implemented within the particlefilter module:

  • selector_deterministic_linear (default): Iteratively selects the sample that has the lowest function value. After selection of a sample, the function value of that point is artificially set to the 2nd highest unique function value.

  • selector_stochastic_uniform: Samples are selected randomly with replacement.

  • selector_stochastic_linear: Samples are selected with replacement, but have different probabilities to be selected. The lowest point has the highest probability to be selected, the highest point is never selected. The function values between this minimum and maximum determine linearly the probabilities for all other points. The unnormalised probability for this selection can be written as:

    q = 1 - (function_value - minimum) / (maximum - minimum)

  • selector_stochastic_softmax: Same as selector_stochastic_linear, but with a softmax probability function.

Survival rate

A part of the samples and function values from each iteration is carried over to each next iteration. The precise points correspond to the best x% of the points (i.e. the x% points with their function values closest to the minimum). The value of x here is configurable via the survival_rate parameter at construction. The default for this parameter is 0.2.

Callback

The ParticleFilter implements a possibility to provide a custom callback function that will be called at each iteration except the first one. This function has to be provided as the callback argument at construction and should have the following signature:

def callback(alg, ind, samples, values, stdevs, width, x, y)

The arguments of this function are:

  • alg: the ParticleFilter object itself. Sounds weird to pass this, but it gives you direct access to all internal properties of the object, together with all object methods.
  • ind: Indices of the samples (and values) of the previous iteration that were used to construct the normal distributions to sample the current iteration from.
  • samples: The samples from the previous iteration. Can be used in combination with the ind argument to find the samples that were taken as the means for the gaussians to sample the current iteration from
  • values: Same as samples, but now the function values for the previous iteration.
  • stdevs: The standard deviations for the gaussians used to sample the current iteration. The stdevs correspond to the means samples[ind].
  • width: The width parameter used in this iteration
  • x: Samples SAMPLED in this iteration. This does not include the samples that were carried over from the previous iteration. These can be collected after the run is over from the samples.csv file.
  • y: Function values of the samples taken this generation as passed through the x argument.

Misc

The class automatically keeps track of the iteration number in the iteration property. This information can be used in for example the width scheduler.

Running the algorithm: example

The following code is an example of how the ParticleFilter can be used. It used the TestFunction defined above.

# Define ingredients
algorithm = ParticleFilter(seed_size=100, iteration_size=100, width_decay=0.8, hard_ranges=True)
function = ExampleFunction()

# Run algorithm on function
algorithm.check_testfunction(function)
x0, y0 = algorithm(function)
for _ in range(25):
    _, _ = algorithm(function)
xn, yn = algorithm(function)
Clone this wiki locally