This example illustrates how to perform Kriging regression on noisy data. The example is based on a simple one-dimensional function similar to the one in uq_Example_Kriging_01_1D
, now with additive Gaussian noise added to the model response. Three cases of noise variance are considered: unknown homogeneous (homoscedastic), known homoscedastic, and known heterogeneous (heteroscedastic).
from uqpylab import sessions, display_util
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
# Initialize common plotting parameters
display_util.load_plt_defaults()
uq_colors = display_util.get_uq_color_order()
# Start the session
mySession = sessions.cloud()
# (Optional) Get a convenient handle to the command line interface
uq = mySession.cli
# Reset the session
mySession.reset()
uqpylab.sessions :: INFO :: A new session (8bd3e32439a144f7a56210610c081eb1) started. uqpylab.sessions :: INFO :: Reset successful.
uq.rng(0,'twister');
np.random.seed(0)
The computational model is a simple analytical function defined by:
$$y(x) = x \sin(x), \; x \in [-3\pi, 3\pi]$$Specify this model using a string and create a MODEL object:
ModelOpts = {
'Type': 'Model',
'ModelFun': 'branin.model',
}
myModel = uq.createModel(ModelOpts)
The probabilistic input model consists of two independent uniform random variables:
$$X_1 \in \mathcal{U}(-5, 10), \, X_2 \in \mathcal{U}(0, 15)$$Specify the marginals and create a UQCloud INPUT object:
InputOpts = {
'Marginals': [{
'Type': 'Uniform',
'Parameters': [-5, 10]
},
{
'Type': 'Uniform',
'Parameters': [0, 15]
}
]
}
myInput = uq.createInput(InputOpts)
Select the metamodeling tool and the Kriging module:
MetaOpts = {
'Type': 'Metamodel',
'MetaType': 'Kriging'
}
Assign the previously created INPUT and MODEL objects to the metamodel options:
MetaOpts['Input'] = myInput['Name']
MetaOpts['FullModel'] = myModel['Name']
If an INPUT object and a MODEL object are specified, an experimental design is automatically generated and the corresponding model responses are computed.
Specify the sampling strategy and the number of sample points for the experimental design:
MetaOpts['ExpDesign'] = {
'Sampling': 'LHS',
'NSamples': 35
}
Use a smooth correlation family, e.g., Matern 5/2:
MetaOpts['Corr'] = {
'Family': 'Matern-5_2'
}
Create the Kriging metamodel:
myKriging = uq.createModel(MetaOpts)
uqpylab.sessions :: INFO :: Received intermediate compute request, function: branin.model. uqpylab.sessions :: INFO :: Carrying out local computation... uqpylab.sessions :: INFO :: Local computation complete. uqpylab.sessions :: INFO :: Starting transmission of intermediate compute results ((35, 1))... uqpylab.sessions :: INFO :: Intermediate compute results sent.
Print out a report on the resulting Kriging object:
uq.print(myKriging)
%---------------- Kriging metamodel ----------------% Object Name: Model 2 Input Dimension: 2 Output Dimension: 1 Experimental Design X size: [35x2] Y size: [35x1] Sampling: LHS Trend Type: ordinary Degree: 0 Beta: [ 517.09721 ] Gaussian Process (GP) Corr. type: ellipsoidal Corr. isotropy: anisotropic Corr. family: Matern-5_2 sigma^2: 3.73821e+04 Estimation method: Cross-validation (CV) Hyperparameters theta: [ 2.85920 9.99980 ] Optim. method: Hybrid Genetic Alg. GP Regression Mode: interpolation Error estimates Leave-one-out: 7.67774e-04 %---------------------------------------------------%
Create a validation sample of size $10^3$:
Xval = uq.getSample(myInput,1e3)
Evaluate the full model responses at the validation set points:
Yval = uq.evalModel(myModel, Xval)
Evaluate the Kriging predictor mean at the validation set points:
YKRGmean = uq.evalModel(myKriging, Xval)
Visualize the full model responses:
x1, x2 = np.meshgrid(np.linspace(-5,10,20), np.linspace(0,15,20))
Xplot = np.hstack((x1.reshape(-1,1),x2.reshape(-1,1)))
Yplot = uq.evalModel(myModel, Xplot)
[muYKRGplot,varYKRGplot] = uq.evalModel(myKriging,Xplot,nargout=2)
XED = np.array(myKriging['ExpDesign']['X'])
norm = matplotlib.colors.Normalize(vmin=0,vmax=315)
f, axs = plt.subplots(2,2,figsize=(9,9))
im = axs[0][0].contourf(Xplot[:,0].reshape(20,20),Xplot[:,1].reshape(20,20),Yplot.reshape(20,20), levels=20, cmap='Purples',norm=norm)
axs[0][0].set_title(r'$\mathrm{\mathcal{M}({\bf X})}$')
axs[0][0].set_xlabel(r'$\mathrm{{X_1}}$')
axs[0][0].set_ylabel(r'$\mathrm{{X_2}}$')
f.colorbar(im, ax=axs[0][0])
im = axs[0][1].contourf(Xplot[:,0].reshape(20,20),Xplot[:,1].reshape(20,20),muYKRGplot.reshape(20,20), levels= 20, cmap='Purples', norm=norm)
axs[0][1].autoscale(False)
axs[0][1].set_title(r'$\mathrm{\mu_{\widehat{Y}}}^{KRG}$')
axs[0][1].set_xlabel(r'$\mathrm{{X_1}}$')
axs[0][1].set_ylabel(r'$\mathrm{{X_2}}$')
f.colorbar(im, ax=axs[0][1])
axs[0][1].scatter(XED[:,0],XED[:,1],color=uq_colors[0],s=3)
axs[1][0].axis('off')
im = axs[1][1].contourf(Xplot[:,0].reshape(20,20),Xplot[:,1].reshape(20,20),varYKRGplot.reshape(20,20), levels=20, cmap='Purples')
axs[1][1].autoscale(False)
axs[1][1].set_title(r'$\mathrm{\sigma^2_{\widehat{Y}}}^{KRG}$')
axs[1][1].set_xlabel(r'$\mathrm{{X_1}}$')
axs[1][1].set_ylabel(r'$\mathrm{{X_2}}$')
f.colorbar(im, ax=axs[1][1])
axs[1][1].scatter(XED[:,0],XED[:,1],color=uq_colors[0],s=3);
plt.show()
To visually assess the performance of the Kriging metamodel, create a scatter plot of the Kriging predictions (i.e., the mean) vs. the true responses on the validation set:
plt.scatter(Yval, YKRGmean, marker='o',alpha=.7, s=5)
plt.plot([np.min(Yval), np.max(Yval)], [np.min(Yval), np.max(Yval)], 'k')
plt.grid('True');
plt.xlabel('$\mathrm{Y^{true}}$');
plt.ylabel('$\mathrm{\mu_{\widehat{Y}}}^{KRG}$');
plt.tick_params(axis='both')
plt.xlim(np.min(Yval), np.max(Yval))
plt.ylim(np.min(Yval), np.max(Yval))
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
mySession.quit()