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