KRIGING METAMODELING: VARIOUS METHODS¶

This example showcases how to perform Kriging metamodeling for a simple one-dimensional function, using various hyperparameter estimation and optimization methods.

Package imports and initialization¶

In [1]:
from uqpylab import sessions, display_util
import numpy as np
import matplotlib.pyplot as plt

# Initialize common plotting parameters
display_util.load_plt_defaults()
uq_colors = display_util.get_uq_color_order()

Start a remote UQCloud session¶

In [2]:
# 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 (4c3ccb4810574e4ab022ef18267c50cb) started.
 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

In [3]:
uq.rng(0,'twister');

Computational model¶

The computational model is a simple analytical function defined by:

$$ y = x \sin(x), \; x \in [0,15] $$

Specify this model in UQCloud using a string with vectorized operation:

In [4]:
ModelOpts = {
    'Type': 'Model',
    'Name': 'XsinX',
    'ModelFun': 'xsinx.model'
}

myModel = uq.createModel(ModelOpts)

Probabilistic input model¶

The probabilistic input model consists of one uniform random variable: $X \sim \mathcal{U}(0, 15)$

Specify the marginal and create an INPUT object:

In [6]:
InputOpts = {
    'Marginals': {
        'Type': 'Uniform',
        'Parameters': [0, 15]
    }
}

myInput = uq.createInput(InputOpts)

Experimental design and model responses¶

Generate $10$ sample points of $X$ using the Sobol' sequence sampling:

In [7]:
X = uq.getSample(myInput, 10, 'Sobol')

Evaluate the corresponding model responses:

In [8]:
Y = uq.evalModel(myModel,X)

Kriging Metamodels¶

Select the metamodeling tool and the Kriging module:

In [9]:
MetaOpts = {
    'Type': 'Metamodel',
    'MetaType': 'Kriging'
}

Assign the experimental design and the corresponding model responses:

In [10]:
MetaOpts['ExpDesign'] = {
    'X': X.tolist(),
    'Y': Y.tolist()
}

Perform ordinary Kriging (i.e., Kriging with a constant but unknown trend):

In [11]:
MetaOpts['Trend'] = {
    'Type': 'ordinary'
}

Set the correlation family to Matérn 5/2:

In [12]:
MetaOpts['Corr'] = {
    'Family': 'matern-5_2'
}

The previous options are fixed, i.e., all the Kriging metamodels that will be created will use them. The selected optimization and hyperparameter estimation method, however, varies for each Krigin metamodel that is created next:

  • Estimation method: Maximum-Likelihood, Optimization method: BFGS
In [13]:
MetaOpts['EstimMethod'] = 'ML'
MetaOpts['Optim'] = {
    'Method': 'BFGS'
}
myKriging_ML_BFGS = uq.createModel(MetaOpts)
  • Estimation method: Maximum-Likelihood, Optimization method: HGA
In [14]:
MetaOpts['Optim'] = {
    'Method': 'HGA'
}
myKriging_ML_HGA = uq.createModel(MetaOpts)
  • Estimation method: Cross-Validation, Optimization method: BFGS
In [15]:
MetaOpts['EstimMethod'] = 'CV'
MetaOpts['Optim'] = {
    'Method': 'BFGS'
}
myKriging_CV_BFGS = uq.createModel(MetaOpts)
  • Estimation method: Cross-Validation, Optimization method: HGA
In [16]:
MetaOpts['EstimMethod'] = 'CV'
MetaOpts['Optim'] = {
    'Method': 'HGA'
}
myKriging_CV_HGA = uq.createModel(MetaOpts)

Comparison of the metamodels¶

Create a validation set of size $10^3$ over a regular grid:

In [17]:
Xval = uq.getSample(myInput, 10**3, 'grid')

Evaluate the true model response at the validation sample points:

In [18]:
Yval = uq.evalModel(myModel,Xval)

Evaluate the Kriging surrogate predictions on the validation set. In each case, both the mean and the variance of the Kriging predictor are calculated:

In [19]:
Ymu_ML_BFGS, Yvar_ML_BFGS = uq.evalModel(myKriging_ML_BFGS, Xval, nargout=2)
Ymu_ML_HGA, Yvar_ML_HGA = uq.evalModel(myKriging_ML_HGA, Xval, nargout=2)
Ymu_CV_BFGS, Yvar_CV_BFGS = uq.evalModel(myKriging_CV_BFGS, Xval, nargout=2)
Ymu_CV_HGA, Yvar_CV_HGA = uq.evalModel(myKriging_CV_HGA, Xval, nargout=2)

Lastly, comparison plots for the mean and the variance of the Kriging predictors are generated. They are divided into two groups based on the hyperparameter estimation method.

  • Maximum likelihood estimation (mean):
In [25]:
plt.plot(Xval, Yval, '-k',linewidth=1.0)
plt.plot(Xval, Ymu_ML_BFGS,'-',c=uq_colors[1])
plt.plot(Xval, Ymu_ML_HGA, '--',c=uq_colors[0])
plt.plot(X, Y, 'ko',markersize=2)

plt.xlim([0, 15])
plt.ylim([-15, 20])
plt.legend(['True model', 'Kriging, optim.method: BFGS',
           'Kriging, optim. method: HGA', 'Observations'],
          loc='upper left')
plt.xlabel('$\mathrm{X}$')
plt.ylabel('$\mathrm{Y}$')
plt.title('Maximum likelihood - Mean');
plt.show()
  • Maximum likelihood estimation (variance):
In [26]:
plt.plot(X, np.zeros(X.shape), 'k', linewidth=1)
plt.plot(Xval, Yvar_ML_BFGS, '-', c=uq_colors[1])
plt.plot(Xval, Yvar_ML_HGA, '--', c=uq_colors[0])
plt.plot(X, np.zeros(X.shape), 'ko', markersize=2)
plt.xlim([0, 15])
plt.ylim([0, 25])
plt.legend(['True model', 'Kriging, optim.method: BFGS',
           'Kriging, optim. method: HGA', 'Observations'],
          loc='upper center')
plt.xlabel('$\mathrm{X}$')
plt.ylabel('$\mathrm{\sigma^2_{\widehat{Y}}}$')
plt.title('Maximum likelihood - Variance');
plt.show()
  • Cross-validation-based estimation (mean):
In [30]:
plt.plot(Xval, Yval, 'k', linewidth=1)
plt.plot(Xval, Ymu_CV_BFGS, '-',c=uq_colors[1])
plt.plot(Xval, Ymu_CV_HGA, '--',c=uq_colors[0])
plt.plot(X, Y, 'ko', markersize=3)

plt.xlim([0, 15])
plt.ylim([-15, 25])
plt.legend(['True model', 'Kriging, optim.method: BFGS',
           'Kriging, optim. method: HGA', 'Observations'],
          loc='upper center')
plt.xlabel('$\mathrm{X}$')
plt.ylabel('$\mathrm{Y}$')
plt.title('Cross-validation - Mean');
plt.show()
  • Cross-validation-based estimation (variance):
In [29]:
plt.plot(Xval, np.zeros(Yval.shape), 'k', linewidth=1)
plt.plot(Xval, Yvar_CV_BFGS, '-.',c=uq_colors[1])
plt.plot(Xval, Yvar_CV_HGA, '--',c=uq_colors[0])
plt.plot(X, np.zeros(X.shape), 'ko')

plt.xlim([0, 15])
plt.ylim([0, 25])

plt.legend(['True model', 'Kriging, optim.method: BFGS',
           'Kriging, optim. method: HGA', 'Observations'],
          loc='upper center')
plt.xlabel('$\mathrm{X}$')
plt.ylabel('$\mathrm{\sigma^2_{\widehat{Y}}}$')
plt.title('Cross-validation - Variance');
plt.show()

Terminate the remote UQCloud session¶

In [ ]:
mySession.quit()