# Regularisation#

We can maintain particle diversity with the **post-regularised particle filter**, where resampling draws particles from a continuous distribution, rather than directly from the particle ensemble.
The continuous distribution is constructed using a multi-variate Gaussian kernel, whose covariance matrix is the normalised particle covariance.
This means that when the particles are resampled, we will obtain many particles that are **similar, but not identical** to the particles with the largest weights (i.e., those that best describe the past observations).

As demonstrated in the figure below, the post-regularised particle filter greatly improves the model fit and forecasts. However, due to the chaotic nature of the system, the forecasts for \(x(t)\) and \(y(t)\) become quite uncertain by \(t=25\).

In order to use the post-regularised particle filter, we need to:

Identify which fields in the particle state vector can be smoothed by the filter. This is defined by the simulation model’s

`can_smooth()`

method.Enable regularisation by setting

`filter.regularisation.enabled = true`

in the scenario definition.Define the lower and upper bounds for each field that should be smoothed (

`filter.regularisation.bounds.FIELD`

) in the scenario definition.

Note

The post-regularised particle filter will only smooth fields that are returned by the `can_smooth()`

method **and** for which lower and/or upper bounds have been defined.

You can omit the lower or upper bounds, if they are not appropriate.
To smooth a field without imposing any bounds, define an empty table for that field.
See the bounds for `z`

in the scenario definition below.

```
class Lorenz63(OdeModel):
def field_types(self, ctx):
r"""
Define the state vector :math:`[\sigma, \rho, \beta, x, y, z]^T`.
"""
return [
('sigma', float),
('rho', float),
('beta', float),
('x', float),
('y', float),
('z', float),
]
def d_dt(self, time, xt, ctx, is_forecast):
"""
The right-hand side of the ODE system.
:param time: The current time.
:param xt: The particle state vectors.
:param ctx: The simulation context.
:param is_forecast: True if this is a forecasting simulation.
"""
rates = np.zeros(xt.shape, xt.dtype)
rates['x'] = xt['sigma'] * (xt['y'] - xt['x'])
rates['y'] = xt['x'] * (xt['rho'] - xt['z']) - xt['y']
rates['z'] = xt['x'] * xt['y'] - xt['beta'] * xt['z']
return rates
def can_smooth(self):
"""Indicate which state vector fields can be smoothed."""
return {'sigma', 'rho', 'beta', 'x', 'y', 'z'}
```

```
1# NOTE: Save this file as 'lorenz63_forecast_regularised.toml'
2
3[components]
4model = "pypfilt.examples.lorenz.Lorenz63"
5time = "pypfilt.Scalar"
6sampler = "pypfilt.sampler.LatinHypercube"
7summary = "pypfilt.summary.HDF5"
8
9[time]
10start = 0.0
11until = 25.0
12steps_per_unit = 10
13summaries_per_unit = 10
14
15[prior]
16sigma = { name = "constant", args.value = 10 }
17rho = { name = "constant", args.value = 28 }
18beta = { name = "constant", args.value = 2.66667 }
19x = { name = "uniform", args.loc = -5, args.scale = 10 }
20y = { name = "uniform", args.loc = -5, args.scale = 10 }
21z = { name = "uniform", args.loc = -5, args.scale = 10 }
22
23[summary.tables]
24forecasts.component = "pypfilt.summary.PredictiveCIs"
25forecasts.credible_intervals = [50, 60, 70, 80, 90, 95]
26sim_z.component = "pypfilt.summary.SimulatedObs"
27sim_z.observation_unit = "z"
28
29[observations.x]
30model = "pypfilt.examples.lorenz.ObsLorenz63"
31file = "lorenz63-x.ssv"
32
33[observations.y]
34model = "pypfilt.examples.lorenz.ObsLorenz63"
35file = "lorenz63-y.ssv"
36
37[observations.z]
38model = "pypfilt.examples.lorenz.ObsLorenz63"
39file = "lorenz63-z.ssv"
40
41[filter]
42particles = 500
43prng_seed = 2001
44history_window = -1
45resample.threshold = 0.25
46regularisation.enabled = true
47
48[filter.regularisation.bounds]
49x = { min = -50, max = 50 }
50y = { min = -50, max = 50 }
51z = {}
52
53[scenario.forecast_regularised]
```

Note

Call `save_lorenz63_scenario_files()`

to save this scenario file (and the others used in this tutorial) in the working directory.

```
def run_lorenz63_forecast_regularised(filename=None):
scenario_file = 'lorenz63_forecast_regularised.toml'
instances = list(pypfilt.load_instances(scenario_file))
instance = instances[0]
# Run a forecast from t = 20.
forecast_time = 20
context = instance.build_context()
return pypfilt.forecast(context, [forecast_time], filename=filename)
```