Example 04: noisy sequential prediction

Script: examples/example04_sequential_prediction_with_noise.py

Purpose

The script shows how to use SequentialPrediction with noisy scalar observations. The model class is Model_Noisy_ConstantMean_Maternp_REML. The physical input is augmented with one observation-noise variance column, so observations and prediction points carry different noise values. This is the standard noisy-observation form of GP prediction, with the noise variance added to the observation covariance [2, 7].

What is computed

• noisy observations of the one-dimensional twobumps function.

• covariance parameters selected by noisy REML.

• posterior means and variances at zero-noise prediction points.

• sequential additions selected by maximum posterior variance.

• conditional sample paths from the final noisy model.

Main objects

• gpmpcontrib.computerexperiment.StochasticComputerExperiment

• gpmpcontrib.Model_Noisy_ConstantMean_Maternp_REML

• gpmpcontrib.SequentialPrediction

• gpmpcontrib.plot.plot_1d

Outputs

Run python examples/example04_sequential_prediction_with_noise.py from the repository root to execute the example. Regenerate the static figure with cd docs && python make_example_results.py.

Red points are noisy observations. Several observations are replicated at the same physical input. The prediction grid is passed with zero observation noise, so the posterior variance represents uncertainty on the latent function, not observation noise at a new noisy measurement.

Source excerpt

# -- 1. Define a problem --
input_dim = 1
input_box = [[-1], [1]]
output_dim = 1
noise_variance = 0.2**2

problem = gpc.computerexperiment.StochasticComputerExperiment(
    input_dim,
    input_box,
    single_function=gp.misc.testfunctions.twobumps,
    simulated_noise_variance=noise_variance,
)

# -- 2. Create initial dataset --
nt = 2000
xt = gp.misc.designs.regulargrid(problem.input_dim, nt, problem.input_box)
zt = problem(xt, simulate_noise=False)

ind = [100, 1600] + [1000] * 15
xi = xt[ind]
zi = problem(xi)

# -- 3. Create a Model and a SequentialPrediction object --
model = gpc.Model_Noisy_ConstantMean_Maternp_REML(
    "Simple function",
    output_dim=problem.output_dim,
    mean_specification={"type": "constant"},
    covariance_specification={"p": 3},
)

# Using the append_noise_variance function to manage noise data
xi_with_noise_variance = append_noise_variance(xi, noise_variance, problem.output_dim)
xt_with_zero_noise_variance = append_noise_variance(xt, 0, problem.output_dim)

sp = gpc.SequentialPrediction(model)
sp.set_data_with_model_selection(xi_with_noise_variance, zi)

zpm, zpv = sp.predict(xt_with_zero_noise_variance)

title = f"1D GP model with {sp.ni} observations"
plot_1d(xt, zt, xi, zi, zpm, zpv, title=title)