KRIGING METAMODELING: TREND TYPES

This example showcases how to perform Kriging for a simple one-dimensional function using various trend types.

Package imports

from uqpylab import sessions, display_util
import numpy as np
import matplotlib.pyplot as plt

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

Start a remote UQCloud session

# 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 (2c52093ebe194e49b2206b35444470dd) started.
 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility

uq.rng(100,'twister');

Computational model

The computational model is a simple analytical function defined by:

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

Specify this model using a string and create a MODEL object:

ModelOpts = {
    'Type': 'Model',
    'ModelFun': 'xsinx.model',
    'isVectorized': 'true'
}

myModel = uq.createModel(ModelOpts)

Probabilistic input model

The probabilistic input model consists of a single uniform random variable:

$X \sim \mathcal{U}(0, 15)$

Specify its marginal distribution and create a INPUT object:

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:

X = uq.getSample(N=10, Method='Sobol')

Y = uq.evalModel(myModel, X)

Kriging metamodels

Three trend types are considered in this example:

• Constant (so-called ordinary Kriging)
• 3rd-degree polynomial
• Custom type

First, specify all the options common among all the three Kriging metamodels.

Select the metamodeling tool and the Kriging module:

MetaOpts = {
    "Type": "Metamodel",
    "MetaType": "Kriging"
}

Assign the experimental design and the corresponding model responses:

MetaOpts["ExpDesign"] = {
    "X": X.tolist(),
    "Y": Y.tolist()
}

Set the correlaton family to Matérn 3/2:

MetaOpts['Corr'] = {
    "Family": "Matern-3_2"
}

Use cross-validation to estimate the hyperparameters:

MetaOpts['EstimMethod'] = 'CV'

Use a hybrid genetic algorithm (GA)-based optimization to solve the cross-validation problem:

MetaOpts['Optim'] = {
    'Method': 'HGA'
}

In the hybrid genetic algorithm, the final solution of the genetic algorithm is used as the starting point of a gradient-based optimization method (BFGS).

Specify the maximum number of generations, convergence tolerance, and population size of the GA-based optimization:

MetaOpts['Optim']['MaxIter'] = 100
MetaOpts['Optim']['Tol'] = 1e-6
MetaOpts['Optim']['HGA'] = {
    'nPop': 40
}

Ordinary Kriging

Set the Kriging trend type to ordinary Kriging:

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

Create the first Kriging metamodel:

myKrigingOrdinary = uq.createModel(MetaOpts)

Polynomial Trend

A second Kriging metamodel is created using a 3rd-degree polynomial as the trend function, while the other options remain the same.

Set the Kriging trend type to polynomial and specify the degree:

MetaOpts['Trend'] = {
    'Type': 'polynomial',
    'Degree': 3
}

Create the second Kriging metamodel:

myKrigingPolynomial = uq.createModel(MetaOpts)

Custom Trend

A third Kriging metamodel is created with a custom functional basis as the trend:

$$f(x) = x^2 + \sqrt{|x|}$$

Set the Kriging trend type to custom and specify the custom function using a function handle:

MetaOpts['Trend'] = {
    'Type': 'custom',
    'CustomF': '@(X) X.^2 + sqrt(abs(X))'
}

Create the third Kriging metamodel:

myKrigingCustom = uq.createModel(MetaOpts)

VALIDATION

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

Xval = uq.getSample(N=1e3, Method='grid')

Evaluate the true model responses on the validation set:

Yval = uq.evalModel(myModel, Xval)

Evaluate the corresponding responses of each of the three Kriging metamodels:

[YMeanOrdinary,YVarOrdinary] = uq.evalModel(myKrigingOrdinary, Xval, nargout=2)
[YMeanPolynomial,YVarPolynomial] = uq.evalModel(myKrigingPolynomial, Xval, nargout=2)
[YMeanCustom,YVarCustom] = uq.evalModel(myKrigingCustom, Xval, nargout=2)

For each metamodel, the mean and the variance of the Kriging predictor are calculated.

Compare the mean prediction of each Kriging metamodel on the validation set (also taking into account the true model responses):

plt.plot(Xval, Yval, '-k',linewidth=1.0)
plt.plot(Xval, YMeanOrdinary,'-',c=uq_colors[1])
plt.plot(Xval, YMeanPolynomial, '--',c=uq_colors[0])
plt.plot(Xval, YMeanCustom, '--',c='green')
plt.plot(X, Y, 'ko',markersize=2)

plt.xlim([0, 15])
plt.ylim([-15, 25])
plt.legend(['True model', 'Kriging, R: Ordinary Kriging',
           'Kriging, R: 3rd deg. polynomial trend',
           'Kriging, R: Custom trend',
           'Observations'],
          loc='upper left')
plt.xlabel('$\mathrm{X}$')
plt.ylabel('$\mathrm{Y}$')
plt.title('Maximum likelihood - Mean');
plt.show()

Finally, compare the variance that is predicted by each Kriging metamodel:

plt.plot(Xval, YVarOrdinary,'-',c=uq_colors[1])
plt.plot(Xval, YVarPolynomial, '--',c=uq_colors[0])
plt.plot(Xval, YVarCustom, '--',c='green')
plt.plot(X, np.zeros_like(X), 'ko',markersize=2)

plt.xlim([0, 15])
plt.ylim([0, 45])
plt.legend(['Kriging, R: Ordinary Kriging',
           'Kriging, R: 3rd deg. polynomial trend',
           'Kriging, R: Custom trend',
           'Observations'],
          loc='upper left')
plt.xlabel('$\mathrm{X}$')
plt.ylabel('$\mathrm{Y}$')
plt.title('Maximum likelihood - Variance');
plt.show()

Terminate the remote UQCloud session

mySession.quit()