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

```
def run_lorenz63_forecast_regularised():
scenario_file = 'lorenz63_forecast_regularised.toml'
with open(scenario_file, 'w') as f:
f.write(pypfilt.examples.lorenz.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=None)
```