` example.
Definition of the model
-----------------------
.. code-block:: default
import openturns as ot
ot.RandomGenerator.SetSeed(0)
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)
We load the cantilever beam use case
.. code-block:: default
from openturns.usecases import cantilever_beam as cantilever_beam
cb = cantilever_beam.CantileverBeam()
We load the function (model) which evaluates the output Y depending on the inputs.
.. code-block:: default
model = cb.model
Then we define the distribution of the input random vector.
.. code-block:: default
dimension = cb.dim # number of inputs
myDistribution = cb.distribution
Create the design of experiments
--------------------------------
We consider a simple Monte-Carlo sampling as a design of experiments. This is why we generate an input sample using the `getSample` method of the distribution. Then we evaluate the output using the `model` function.
.. code-block:: default
sampleSize_train = 20
X_train = myDistribution.getSample(sampleSize_train)
Y_train = model(X_train)
Create the Legendre basis
-------------------------
We first create a Legendre basis of univariate polynomials. In order to convert then into multivariate polynomials, we use a linear enumerate function.
The `LegendreFactory` class creates Legendre polynomials.
.. code-block:: default
univariateFactory = ot.LegendreFactory()
This factory corresponds to the `Uniform` distribution in the [-1,1] interval.
.. code-block:: default
univariateFactory.getMeasure()
.. raw:: html
Uniform(a = -1, b = 1)
This interval does not correspond to the interval on which the input marginals are defined (we will come back to this topic later), but this will, anyway, create a consistent trend for the kriging.
.. code-block:: default
polyColl = [univariateFactory]*dimension
.. code-block:: default
enumerateFunction = ot.LinearEnumerateFunction(dimension)
productBasis = ot.OrthogonalProductPolynomialFactory(polyColl, enumerateFunction)
.. code-block:: default
functions = []
numberOfTrendCoefficients = 12
for i in range(numberOfTrendCoefficients):
multivariatepolynomial = productBasis.build(i)
print(multivariatepolynomial)
functions.append(multivariatepolynomial)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
1
1.73205 * x0
1.73205 * x1
1.73205 * x2
1.73205 * x3
-1.11803 + 3.3541 * x0^2
(1.73205 * x0) * (1.73205 * x1)
(1.73205 * x0) * (1.73205 * x2)
(1.73205 * x0) * (1.73205 * x3)
-1.11803 + 3.3541 * x1^2
(1.73205 * x1) * (1.73205 * x2)
(1.73205 * x1) * (1.73205 * x3)
.. code-block:: default
basis = ot.Basis(functions)
Create the metamodel
--------------------
In order to create the kriging metamodel, we first select a constant trend with the `ConstantBasisFactory` class. Then we use a squared exponential covariance model. Finally, we use the `KrigingAlgorithm` class to create the kriging metamodel, taking the training sample, the covariance model and the trend basis as input arguments.
.. code-block:: default
covarianceModel = ot.SquaredExponential([1.]*dimension, [1.0])
.. code-block:: default
algo = ot.KrigingAlgorithm(X_train, Y_train, covarianceModel, basis)
algo.run()
result = algo.getResult()
krigingWithConstantTrend = result.getMetaModel()
The `getTrendCoefficients` method returns the coefficients of the trend.
.. code-block:: default
result.getTrendCoefficients()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
[class=Point name=Unnamed dimension=12 values=[6.80011e-46,8.00088e-35,3.49721e-43,2.99508e-45,1.7335e-52,1.05339e-23,4.11637e-32,3.52382e-34,2.03902e-41,2.02866e-40,1.54017e-42,8.91627e-50]]
We see that the number of coefficients in the trend corresponds to the number of functions in the basis.
.. code-block:: default
result.getCovarianceModel()
.. raw:: html
SquaredExponential(scale=[1,1,1,1], amplitude=[0.0316175])
The `SquaredExponential` model has one amplitude coefficient and 4 scale coefficients. This is because this covariance model is anisotropic : each of the 4 input variables is associated with its own scale coefficient.
Create an orthogonal multivariate polynomial factory
----------------------------------------------------
In order to create a Legendre basis which better corresponds to the input marginals, we could consider the orthogonal basis which would be associated to uniform marginals. To compute the bounds of these uniform distributions, we may consider the 1% and 99% quantiles of each marginal.
There is, however, a simpler way to proceed. We can simply orthogonalize the input marginals and create the orthogonal polynomial basis corresponding to the inputs. This corresponds to the method we would use in the polynomial chaos.
We first create the polynomial basis which corresponds to the inputs.
.. code-block:: default
multivariateBasis = ot.OrthogonalProductPolynomialFactory([cb.E, cb.F, cb.L, cb.I])
Then we create the multivariate basis which has maximum degree equal to 2.
.. code-block:: default
totalDegree = 2
enumerateFunction = multivariateBasis.getEnumerateFunction()
numberOfTrendCoefficients = enumerateFunction.getStrataCumulatedCardinal(totalDegree)
numberOfTrendCoefficients
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
15
.. code-block:: default
functions = []
for i in range(numberOfTrendCoefficients):
multivariatepolynomial = productBasis.build(i)
print(multivariatepolynomial)
functions.append(multivariatepolynomial)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
1
1.73205 * x0
1.73205 * x1
1.73205 * x2
1.73205 * x3
-1.11803 + 3.3541 * x0^2
(1.73205 * x0) * (1.73205 * x1)
(1.73205 * x0) * (1.73205 * x2)
(1.73205 * x0) * (1.73205 * x3)
-1.11803 + 3.3541 * x1^2
(1.73205 * x1) * (1.73205 * x2)
(1.73205 * x1) * (1.73205 * x3)
-1.11803 + 3.3541 * x2^2
(1.73205 * x2) * (1.73205 * x3)
-1.11803 + 3.3541 * x3^2
.. code-block:: default
basis = ot.Basis(functions)
.. code-block:: default
algo = ot.KrigingAlgorithm(X_train, Y_train, covarianceModel, basis)
algo.run()
result = algo.getResult()
krigingWithConstantTrend = result.getMetaModel()
The `getTrendCoefficients` method returns the coefficients of the trend.
.. code-block:: default
result.getTrendCoefficients()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
[class=Point name=Unnamed dimension=15 values=[6.80011e-46,8.00088e-35,3.49721e-43,2.99508e-45,1.7335e-52,1.05339e-23,4.11637e-32,3.52382e-34,2.03902e-41,2.02866e-40,1.54017e-42,8.91627e-50,1.39896e-44,7.63508e-52,-7.60276e-46]]
Conclusion
----------
The trend that we have configured corresponds to the basis that we would have used in a full polynomial chaos computed with least squares.
Other extensions of this work would be:
* to use a Fourier basis with `FourierSeriesFactory`,
* wavelets with `HaarWaveletFactory`,
or any other univariate factory.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 0.009 seconds)
.. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_plot_kriging_beam_arbitrary_trend.py:
.. only :: html
.. container:: sphx-glr-footer
:class: sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_kriging_beam_arbitrary_trend.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_kriging_beam_arbitrary_trend.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_