Gaussian Process Regression : cantilever beam model

In this example, we create a Gaussian Process Regression (GPR) metamodel of the cantilever beam. We use a squared exponential covariance kernel for the Gaussian process. In order to estimate the hyper-parameters, we use a design of experiments of size 20.

from openturns.usecases import cantilever_beam
import openturns as ot
import openturns.experimental as otexp
import openturns.viewer as otv

# sphinx_gallery_thumbnail_number = 3

Definition of the model

We load the cantilever beam use case :

cb = cantilever_beam.CantileverBeam()

We define the function which evaluates the output depending on the inputs.

model = cb.model

Then we define the distribution of the input random vector.

myDistribution = cb.distribution

Create the design of experiments

We consider a simple Monte Carlo sample as a design of experiments. This is why we generate an input sample using the method getSample() of the distribution. Then we evaluate the output using the model function.

sampleSize_train = 20
X_train = myDistribution.getSample(sampleSize_train)
Y_train = model(X_train)

The following figure presents the distribution of the vertical deviations Y on the training sample. We observe that the large deviations occur less often.

histo = ot.HistogramFactory().build(Y_train).drawPDF()
histo.setXTitle("Vertical deviation (cm)")
histo.setTitle("Distribution of the vertical deviation")
histo.setLegends([""])
view = otv.View(histo)

Create the metamodel

In order to create the GPR metamodel, we first select a constant trend with the ConstantBasisFactory class. Then we use a squared exponential covariance kernel. The SquaredExponential kernel has one amplitude coefficient and 4 scale coefficients. This is because this covariance kernel is anisotropic : each of the 4 input variables is associated with its own scale coefficient.

basis = ot.ConstantBasisFactory(cb.dim).build()
covarianceModel = ot.SquaredExponential(cb.dim)

Finally, we use the GaussianProcessFitter and GaussianProcessRegression classes to create the GPR metamodel. It requires a training sample, a covariance kernel and a trend basis as input arguments.

fitter_algo = otexp.GaussianProcessFitter(X_train, Y_train, covarianceModel, basis)

The method run() of the class GaussianProcessFitter optimizes the Gaussian process hyperparameters and the method run() of the class GaussianProcessRegression conditions the Gaussian process to the data set.

We can then print the constant trend of the metamodel, estimated using the least squares method.

fitter_algo.run()
fitter_result = fitter_algo.getResult()
gpr_algo = otexp.GaussianProcessRegression(fitter_result)
gpr_algo.run()
gpr_result = gpr_algo.getResult()
gprMetamodel = gpr_result.getMetaModel()

The method getTrendCoefficients() of the class GaussianProcessRegressionResult returns the coefficients of the trend.

print(gpr_result.getTrendCoefficients())

[0.175118]

We can also print the hyperparameters of the covariance model, which have been estimated by maximizing the likelihood.

gpr_result.getCovarianceModel()

SquaredExponential(scale=[1.39249e+10,149.298,0.17976,4.88951e-08], amplitude=[0.0397017])

Validate the metamodel

We finally want to validate the GPR metamodel. This is why we generate a validation sample with size 100 and we evaluate the output of the model on this sample.

sampleSize_test = 100
X_test = myDistribution.getSample(sampleSize_test)
Y_test = model(X_test)

The class MetaModelValidation makes the validation easy. To create it, we use the validation samples and the metamodel.

val = ot.MetaModelValidation(Y_test, gprMetamodel(X_test))

The method computeR2Score() computes the R2 score.

R2 = val.computeR2Score()[0]
print(R2)

0.997701717379125

The residuals are the difference between the model and the metamodel.

r = val.getResidualSample()
graph = ot.HistogramFactory().build(r).drawPDF()
graph.setXTitle("Residuals (cm)")
graph.setTitle("Distribution of the residuals")
graph.setLegends([""])
view = otv.View(graph)

We observe that the negative residuals occur with nearly the same frequency of the positive residuals: this is a first sign of good quality. The method drawValidation() compares the observed outputs and the metamodel outputs.

graph = val.drawValidation()
graph.setTitle("R2 = %.2f%%" % (100 * R2))
view = otv.View(graph)

Display all figures

otv.View.ShowAll()