.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here ` to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_meta_modeling_polynomial_chaos_metamodel_plot_chaos_ishigami.py:
Create a polynomial chaos for the Ishigami function: a quick start guide to polynomial chaos
============================================================================================
In this example, we create a polynomial chaos for the :ref:`Ishigami function`. We create a sparse polynomial with maximum total degree equal to 8.
Define the model
----------------
.. code-block:: default
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
import numpy as np
ot.Log.Show(ot.Log.NONE)
ot.RandomGenerator.SetSeed(0)
We load the Ishigami model :
.. code-block:: default
from openturns.usecases import ishigami_function as ishigami_function
im = ishigami_function.IshigamiModel()
The `IshigamiModel` data class contains the input distribution :math:`X=(X_1, X_2, X_3)` in `im.distributionX` and the Ishigami function in `im.model`.
We also have access to the input variable names with
.. code-block:: default
input_names = im.distributionX.getDescription()
Draw the function
-----------------
Create a training sample
.. code-block:: default
N = 1000
inputSample = im.distributionX.getSample(N)
outputSample = im.model(inputSample)
.. code-block:: default
def plotXvsY(sampleX, sampleY, figsize=(15,3)):
import pylab as pl
import openturns.viewer
dimX = sampleX.getDimension()
inputdescr = sampleX.getDescription()
fig = pl.figure(figsize=figsize)
for i in range(dimX):
ax = fig.add_subplot(1, dimX, i+1)
graph = ot.Graph('', inputdescr[i], 'Y', True, '')
cloud = ot.Cloud(sampleX[:,i],sampleY)
graph.add(cloud)
_ = ot.viewer.View(graph, figure=fig, axes=[ax])
return None
plotXvsY(inputSample, outputSample)
.. image:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_ishigami_001.png
:alt: plot chaos ishigami
:class: sphx-glr-single-img
.. code-block:: default
graph = ot.HistogramFactory().build(outputSample).drawPDF()
view = viewer.View(graph)
.. image:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_ishigami_002.png
:alt: y0 PDF
:class: sphx-glr-single-img
We see that the distribution of the output has two modes.
Create the polynomial chaos model
---------------------------------
Create a training sample
.. code-block:: default
N = 100
inputTrain = im.distributionX.getSample(N)
outputTrain = im.model(inputTrain)
Create the chaos.
We could use only the input and output training samples: in this case, the distribution of the input sample is computed by selecting the distribution that has the best fit.
.. code-block:: default
chaosalgo = ot.FunctionalChaosAlgorithm(inputTrain, outputTrain)
Since the input distribution is known in our particular case, we instead create the multivariate basis from the distribution.
.. code-block:: default
multivariateBasis = ot.OrthogonalProductPolynomialFactory([im.X1, im.X2, im.X3])
selectionAlgorithm = ot.LeastSquaresMetaModelSelectionFactory()
projectionStrategy = ot.LeastSquaresStrategy(inputTrain, outputTrain, selectionAlgorithm)
totalDegree = 8
enumfunc = multivariateBasis.getEnumerateFunction()
P = enumfunc.getStrataCumulatedCardinal(totalDegree)
adaptiveStrategy = ot.FixedStrategy(multivariateBasis, P)
chaosalgo = ot.FunctionalChaosAlgorithm(inputTrain, outputTrain, im.distributionX, adaptiveStrategy, projectionStrategy)
.. code-block:: default
chaosalgo.run()
result = chaosalgo.getResult()
metamodel = result.getMetaModel()
In order to validate the metamodel, we generate a test sample.
.. code-block:: default
n_valid = 1000
inputTest = im.distributionX.getSample(n_valid)
outputTest = im.model(inputTest)
val = ot.MetaModelValidation(inputTest, outputTest, metamodel)
Q2 = val.computePredictivityFactor()[0]
Q2
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
0.9994752470145457
The Q2 is very close to 1: the metamodel is excellent.
.. code-block:: default
graph = val.drawValidation()
graph.setTitle("Q2=%.2f%%" % (Q2*100))
view = viewer.View(graph)
.. image:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_ishigami_003.png
:alt: Q2=99.95%
:class: sphx-glr-single-img
The metamodel has a good predictivity, since the points are almost on the first diagonal.
Compute and print Sobol' indices
--------------------------------
.. code-block:: default
chaosSI = ot.FunctionalChaosSobolIndices(result)
print(chaosSI.summary())
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
input dimension: 3
output dimension: 1
basis size: 21
mean: [3.51725]
std-dev: [3.70631]
------------------------------------------------------------
Index | Multi-indice | Part of variance
------------------------------------------------------------
6 | [0,4,0] | 0.273917
1 | [1,0,0] | 0.191725
5 | [1,0,2] | 0.135136
11 | [0,6,0] | 0.134366
4 | [3,0,0] | 0.120877
9 | [3,0,2] | 0.0850971
3 | [0,2,0] | 0.0253318
10 | [1,0,4] | 0.0122786
------------------------------------------------------------
------------------------------------------------------------
Component | Sobol index | Sobol total index
------------------------------------------------------------
0 | 0.315184 | 0.557148
1 | 0.442823 | 0.442894
2 | 4.76385e-07 | 0.241993
------------------------------------------------------------
We notice the a coefficient with marginal degree equal to 6 has a significant impact on the output variance. Hence, we cannot get a satisfactory polynomial chaos with total degree less that 6.
Draw Sobol' indices
.. code-block:: default
dim_input = im.distributionX.getDimension()
first_order = [chaosSI.getSobolIndex(i) for i in range(dim_input)]
total_order = [chaosSI.getSobolTotalIndex(i) for i in range(dim_input)]
input_names = im.model.getInputDescription()
graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(input_names, first_order, total_order)
view = viewer.View(graph)
.. image:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_ishigami_004.png
:alt: Sobol' indices
:class: sphx-glr-single-img
The variable which has the largest impact on the output is, taking interactions into account, X1.
We see that X1 has interactions with other variables, since the first order indice is less than the total order indice.
At first order, X3 has no interactions with other variables since its first order indice is close to zero.
Computing the accuracy
----------------------
The interesting point with the Ishigami function is that the exact Sobol' indices are known. We can use that property in order to compute the absolute error on the Sobol' indices for the polynomial chaos.
The following function computes the exact mean, variance and Sobol' indices for this function.
.. code-block:: default
def ishigamiSA(a,b):
'''Exact sensitivity indices of the Ishigami function for given a and b.'''
var = 1.0/2 + a**2/8 + b*np.pi**4/5 + b**2*np.pi**8/18
S1 = (1.0/2 + b*np.pi**4/5+b**2*np.pi**8/50)/var
S2 = (a**2/8)/var
S3 = 0
S13 = b**2*np.pi**8/2*(1.0/9-1.0/25)/var
exact = {
'expectation' : a/2,
'variance' : var,
'S1' : (1.0/2 + b*np.pi**4/5+b**2*np.pi**8.0/50)/var,
'S2' : (a**2/8)/var,
'S3' : 0,
'S12' : 0,
'S23' : 0,
'S13' : S13,
'S123' : 0,
'ST1' : S1 + S13,
'ST2' : S2,
'ST3' : S3 + S13
}
return exact
.. code-block:: default
a = 7.
b = 0.1
exact = ishigamiSA(a,b)
exact
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
{'expectation': 3.5, 'variance': 13.844587940719254, 'S1': 0.31390519114781146, 'S2': 0.4424111447900409, 'S3': 0, 'S12': 0, 'S23': 0, 'S13': 0.2436836640621477, 'S123': 0, 'ST1': 0.5575888552099592, 'ST2': 0.4424111447900409, 'ST3': 0.2436836640621477}
To make the comparisons simpler, we gather the results into a list.
.. code-block:: default
S_exact = [exact["S1"],exact["S2"],exact["S3"]]
ST_exact = [exact["ST1"],exact["ST2"],exact["ST3"]]
Then we perform a loop over the input dimension and compute the absolute error on the Sobol' indices.
.. code-block:: default
for i in range(im.dim):
absoluteErrorS = abs(first_order[i]-S_exact[i])
absoluteErrorST = abs(total_order[i]-ST_exact[i])
print("X%d, Abs.Err. on S=%.1e, Abs.Err. on ST=%.1e" % (i+1, absoluteErrorS,absoluteErrorST))
plt.show()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
X1, Abs.Err. on S=1.3e-03, Abs.Err. on ST=4.4e-04
X2, Abs.Err. on S=4.1e-04, Abs.Err. on ST=4.8e-04
X3, Abs.Err. on S=4.8e-07, Abs.Err. on ST=1.7e-03
We see that the indices are correctly estimated with a low accuracy even if we have use only 100 function evaluations. This shows the good performance of the polynomial chaos in this case.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 0.482 seconds)
.. _sphx_glr_download_auto_meta_modeling_polynomial_chaos_metamodel_plot_chaos_ishigami.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_chaos_ishigami.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_chaos_ishigami.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_