Running forecasts#

In order to fit the simulation model to observations and generate forecasts, we need to:

  1. Choose prior distributions for \(x(t)\), \(y(t)\), and \(z(t)\);

  2. Define an input file for each observation model; and

  3. Record summary statistics such as predictive credible intervals for each observation model, and simulated observations for \(z(t)\).

These changes are indicated by the highlighted lines in the following scenario definition:

An example scenario for generating forecasts for the Lorenz-63 system.#
 1# NOTE: Save this file as 'lorenz63_forecast.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[observations.x]
24model = "pypfilt.examples.lorenz.ObsLorenz63"
25file = "lorenz63-x.ssv"
26
27[observations.y]
28model = "pypfilt.examples.lorenz.ObsLorenz63"
29file = "lorenz63-y.ssv"
30
31[observations.z]
32model = "pypfilt.examples.lorenz.ObsLorenz63"
33file = "lorenz63-z.ssv"
34
35[summary.tables]
36forecasts.component = "pypfilt.summary.PredictiveCIs"
37forecasts.credible_intervals = [50, 60, 70, 80, 90, 95]
38sim_z.component = "pypfilt.summary.SimulatedObs"
39sim_z.observation_unit = "z"
40
41[filter]
42particles = 500
43prng_seed = 2001
44history_window = -1
45resample.threshold = 0.25
46
47[scenario.forecast]

Note

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

Observations will be read from these files when generating forecasts with pypfilt.forecast():

def run_lorenz63_forecast(filename=None):
    scenario_file = 'lorenz63_forecast.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)

If you pass a filename to pypfilt.forecast(), all of the summary tables for the estimation pass and each forecasting pass will be saved to that file, as HDF5 data sets.

Note

HDF5 is a file format that allows you to store lots of data tables and related metadata in a single file, and to load these data tables as if they were NumPy arrays. All of the summary tables recorded by pypfilt are NumPy structured arrays. You can explore HDF5 files with the h5py package, which makes it easy to load and store data tables.

Edit on Bitbucket