` to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_meta_modeling_polynomial_chaos_metamodel_plot_functional_chaos_advanced_ctors.py:
Advanced polynomial chaos construction
======================================
In this example we are going to expose advanced elements in the construction of a polynomial chaos algorithm:
- construction of the multivariate orthonormal basis,
- truncature strategy of the multivariate orthonormal basis,
- evaluation strategy of the approximation coefficients.
In this example, we consider the following function :math:`\mathbb{R}^4 \rightarrow \mathbb{R}`:
.. math::
g(\mathbf{x}) = 1+x_1 x_2 + 2 x_3^2+x_4^4
for any :math:`x_1,x_2,x_3,x_4\in\mathbb{R}`.
We assume that the inputs have normal, uniform, gamma and beta distributions :
.. math::
X_1 \sim \mathcal{N}(0,1), \qquad X_2 \sim \mathcal{U}(-1,1), \qquad X_3 \sim \mathcal{G}(2.75,1), \qquad X_4 \sim \mathcal{B}(2.5,1,-1,2),
and :math:`X_1`, :math:`X_2`, :math:`X_3` and :math:`X_4` are independent.
Define the model and the input distribution
-------------------------------------------
.. code-block:: default
from __future__ import print_function
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)
.. code-block:: default
model = ot.SymbolicFunction(['x1', 'x2', 'x3', 'x4'], ['1+x1*x2 + 2*x3^2+x4^4'])
Create a distribution of dimension 4.
.. code-block:: default
distribution = ot.ComposedDistribution([ot.Normal(), ot.Uniform(), ot.Gamma(2.75, 1.0), ot.Beta(2.5, 1.0, -1.0, 2.0)])
.. code-block:: default
inputDimension = distribution.getDimension()
inputDimension
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
4
STEP 1: Construction of the multivariate orthonormal basis
----------------------------------------------------------
Create the univariate polynomial family collection which regroups the polynomial families for each direction.
.. code-block:: default
polyColl = ot.PolynomialFamilyCollection(inputDimension)
We could use the Krawtchouk and Charlier families (for discrete distributions).
.. code-block:: default
polyColl[0] = ot.KrawtchoukFactory()
polyColl[1] = ot.CharlierFactory()
We could also use the automatic selection of the polynomial which corresponds to the distribution: this is done with the `StandardDistributionPolynomialFactory` class.
.. code-block:: default
for i in range(inputDimension):
marginal = distribution.getMarginal(i)
polyColl[i] = ot.StandardDistributionPolynomialFactory(marginal)
In our specific case, we use specific polynomial factories.
.. code-block:: default
polyColl[0] = ot.HermiteFactory()
polyColl[1] = ot.LegendreFactory()
polyColl[2] = ot.LaguerreFactory(2.75)
# Parameter for the Jacobi factory : 'Probabilty' encoded with 1
polyColl[3] = ot.JacobiFactory(2.5, 3.5, 1)
Create the enumeration function.
The first possibility is to use the `LinearEnumerateFunction`.
.. code-block:: default
enumerateFunction = ot.LinearEnumerateFunction(inputDimension)
Another possibility is to use the `HyperbolicAnisotropicEnumerateFunction`, which gives less weight to interactions.
.. code-block:: default
q = 0.4
enumerateFunction_1 = ot.HyperbolicAnisotropicEnumerateFunction(inputDimension, q)
Create the multivariate orthonormal basis which is the the cartesian product of the univariate basis.
.. code-block:: default
multivariateBasis = ot.OrthogonalProductPolynomialFactory(polyColl, enumerateFunction)
Ask how many polynomials have total degrees equal to k=5.
.. code-block:: default
k = 5
enumerateFunction.getStrataCardinal(k)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
56
Ask how many polynomials have degrees lower or equal to k=5.
.. code-block:: default
enumerateFunction.getStrataCumulatedCardinal(k)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
126
Give the k-th term of the multivariate basis. To calculate its degree, add the integers.
.. code-block:: default
k = 5
enumerateFunction(k)
.. raw:: html
[2,0,0,0]
Build a term of the basis as a Function. Generally, we do not need to construct manually any term, all terms are constructed automatically by a strategy of construction of the basis.
.. code-block:: default
i = 5
Psi_i = multivariateBasis.build(i)
Psi_i
.. raw:: html
-0.707107 + 0.707107 * x0^2
Get the measure mu associated to the multivariate basis.
.. code-block:: default
distributionMu = multivariateBasis.getMeasure()
distributionMu
.. raw:: html
ComposedDistribution(Normal(mu = 0, sigma = 1), Uniform(a = -1, b = 1), Gamma(k = 3.75, lambda = 1, gamma = 0), Beta(alpha = 2.5, beta = 1, a = -1, b = 1), IndependentCopula(dimension = 4))
STEP 2: Truncature strategy of the multivariate orthonormal basis
-----------------------------------------------------------------
FixedStrategy : all the polynomials af degree lower or equal to 2 which corresponds to the 15 first ones.
.. code-block:: default
p = 15
truncatureBasisStrategy = ot.FixedStrategy(multivariateBasis, p)
SequentialStrategy : among the maximumCardinalBasis = 100 first polynomials of the multivariate basis those verfying the convergence criterion.
.. code-block:: default
maximumCardinalBasis = 100
truncatureBasisStrategy_1 = ot.SequentialStrategy(
multivariateBasis, maximumCardinalBasis)
CleaningStrategy : among the maximumConsideredTerms = 500 first polynomials, those which have the mostSignificant = 50 most significant contributions with significance criterion significanceFactor equal to :math:`10^{-4}`
The `True` boolean indicates if we are interested in the online monitoring of the current basis update (removed or added coefficients).
.. code-block:: default
maximumConsideredTerms = 500
mostSignificant = 50
significanceFactor = 1.0e-4
truncatureBasisStrategy_2 = ot.CleaningStrategy(
multivariateBasis, maximumConsideredTerms, mostSignificant, significanceFactor, True)
STEP 3: Evaluation strategy of the approximation coefficients
-------------------------------------------------------------
The technique illustrated is the Least Squares technique where the points come from an design of experiments. Here : the Monte Carlo sampling technique.
.. code-block:: default
sampleSize = 100
evaluationCoeffStrategy = ot.LeastSquaresStrategy(
ot.MonteCarloExperiment(sampleSize))
You can specify the approximation algorithm. This is the algorithm that generates a sequence of basis using Least Angle Regression.
.. code-block:: default
basisSequenceFactory = ot.LARS()
This algorithm estimates the empirical error on each sub-basis using Leave One Out strategy.
.. code-block:: default
fittingAlgorithm = ot.CorrectedLeaveOneOut()
# Finally the metamodel selection algorithm embbeded in LeastSquaresStrategy
approximationAlgorithm = ot.LeastSquaresMetaModelSelectionFactory(
basisSequenceFactory, fittingAlgorithm)
evaluationCoeffStrategy_2 = ot.LeastSquaresStrategy(
ot.MonteCarloExperiment(sampleSize), approximationAlgorithm)
Try integration.
.. code-block:: default
marginalDegrees = [2] * inputDimension
evaluationCoeffStrategy_3 = ot.IntegrationStrategy(
ot.GaussProductExperiment(distributionMu, marginalDegrees))
STEP 4: Creation of the Functional Chaos Algorithm
--------------------------------------------------
The `FunctionalChaosAlgorithm` class combines
* the model : `model`
* the distribution of the input random vector : `distribution`
* the truncature strategy of the multivariate basis
* and the evaluation strategy of the coefficients
.. code-block:: default
polynomialChaosAlgorithm = ot.FunctionalChaosAlgorithm(
model, distribution, truncatureBasisStrategy, evaluationCoeffStrategy)
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 0.006 seconds)
.. _sphx_glr_download_auto_meta_modeling_polynomial_chaos_metamodel_plot_functional_chaos_advanced_ctors.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_functional_chaos_advanced_ctors.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_functional_chaos_advanced_ctors.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_