Example Usage

This example is a reproduction of Figure 1 from http://doi.org/10.1098/rspa.2007.1900 :

The original figure:

https://royalsocietypublishing.org/cms/asset/efa57e07-5384-4503-8b2b-ccbe632ffe87/3251fig1.jpg

Code to reproduce the above figure as close as possible:

 1# Typical imports: Matplotlib, numpy, sklearn and of course our mf2 package
 2import matplotlib.pyplot as plt
 3import mf2
 4import numpy as np
 5from sklearn.gaussian_process import GaussianProcessRegressor as GPR
 6from sklearn.gaussian_process import kernels
 7
 8# Setting up
 9low_x = np.linspace(0, 1, 11).reshape(-1, 1)
10high_x = low_x[[0, 4, 6, 10]]
11diff_x = high_x
12
13low_y = mf2.forrester.low(low_x)
14high_y = mf2.forrester.high(high_x)
15scale = 1.87  # As reported in the paper
16diff_y = np.array([(mf2.forrester.high(x) - scale * mf2.forrester.low(x))[0]
17                   for x in diff_x])
18
19# Training GP models
20kernel = kernels.ConstantKernel(constant_value=1.0) \
21         * kernels.RBF(length_scale=1.0, length_scale_bounds=(1e-1, 10.0))
22
23gp_direct = GPR(kernel=kernel).fit(high_x, high_y)
24gp_low = GPR(kernel=kernel).fit(low_x, low_y)
25gp_diff = GPR(kernel=kernel).fit(diff_x, diff_y)
26
27# Using a simple function to combine the two models
28def co_y(x):
29    return scale * gp_low.predict(x) + gp_diff.predict(x)
30
31# And finally recreating the plot
32plot_x = np.linspace(start=0, stop=1, num=501).reshape(-1, 1)
33plt.figure(figsize=(6, 5), dpi=600)
34plt.plot(plot_x, mf2.forrester.high(plot_x), linewidth=2, color='black', label='$f_e$')
35plt.plot(plot_x, mf2.forrester.low(plot_x), linewidth=2, color='black', linestyle='--',
36         label='$f_c$')
37plt.scatter(high_x, high_y, marker='o', facecolors='none', color='black', label='$y_e$')
38plt.scatter(low_x, low_y, marker='s', facecolors='none', color='black', label='$y_c$')
39plt.plot(plot_x, gp_direct.predict(plot_x), linewidth=1, color='black', linestyle='--',
40         label='kriging through $y_e$')
41plt.plot(plot_x, co_y(plot_x), linewidth=1, color='black', label='co-kriging')
42plt.xlim([0, 1])
43plt.ylim([-10, 20])
44plt.xlabel('$x$')
45plt.ylabel('$y$')
46plt.legend(loc=2)
47plt.tight_layout()
48plt.savefig('../_static/recreating-forrester-2007.png')
49plt.show()

Reproduced figure:

_images/recreating-forrester-2007.png