.. 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_reliability_sensitivity_sensitivity_analysis_plot_functional_chaos_sensitivity.py:
Sobol' sensitivity indices from chaos
=====================================
In this example we are going to compute global sensitivity indices from a functional chaos decomposition.
We study the Borehole function that models water flow through a borehole:
.. math::
\frac{2 \pi T_u (H_u - H_l)}{\ln{r/r_w}(1+\frac{2 L T_u}{\ln(r/r_w) r^2_w K_w}\frac{T_u}{T_l})}
With parameters:
- :math:`r_w`: radius of borehole (m)
- :math:`r`: radius of influence (m)
- :math:`T_u`: transmissivity of upper aquifer (:math:`m^2/yr`)
- :math:`H_u`: potentiometric head of upper aquifer (m)
- :math:`T_l`: transmissivity of lower aquifer (:math:`m^2/yr`)
- :math:`H_l`: potentiometric head of lower aquifer (m)
- :math:`L`: length of borehole (m)
- :math:`K_w`: hydraulic conductivity of borehole (:math:`m/yr`)
.. code-block:: default
from __future__ import print_function
import openturns as ot
from operator import itemgetter
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)
borehole model
.. code-block:: default
dimension = 8
input_names = ['rw', 'r', 'Tu', 'Hu', 'Tl', 'Hl', 'L', 'Kw']
model = ot.SymbolicFunction(input_names,
['(2*pi_*Tu*(Hu-Hl))/(ln(r/rw)*(1+(2*L*Tu)/(ln(r/rw)*rw^2*Kw)+Tu/Tl))'])
coll = [ot.Normal(0.1, 0.0161812),
ot.LogNormal(7.71, 1.0056),
ot.Uniform(63070.0, 115600.0),
ot.Uniform(990.0, 1110.0),
ot.Uniform(63.1, 116.0),
ot.Uniform(700.0, 820.0),
ot.Uniform(1120.0, 1680.0),
ot.Uniform(9855.0, 12045.0)]
distribution = ot.ComposedDistribution(coll)
distribution.setDescription(input_names)
Freeze r, Tu, Tl from model to go faster
.. code-block:: default
selection = [1,2,4]
complement = ot.Indices(selection).complement(dimension)
distribution = distribution.getMarginal(complement)
model = ot.ParametricFunction(model, selection, distribution.getMarginal(selection).getMean())
input_names_copy = list(input_names)
input_names = itemgetter(*complement)(input_names)
dimension = len(complement)
design of experiment
.. code-block:: default
size = 1000
X = distribution.getSample(size)
Y = model(X)
create a functional chaos model
.. code-block:: default
algo = ot.FunctionalChaosAlgorithm(X, Y)
algo.run()
result = algo.getResult()
print(result.getResiduals())
print(result.getRelativeErrors())
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
[0.00224141]
[8.8431e-09]
Quick summary of sensitivity analysis
.. code-block:: default
sensitivityAnalysis = ot.FunctionalChaosSobolIndices(result)
print(sensitivityAnalysis.summary())
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
input dimension: 5
output dimension: 1
basis size: 44
mean: [74.1358]
std-dev: [28.7844]
------------------------------------------------------------
Index | Multi-indice | Part of variance
------------------------------------------------------------
1 | [1,0,0,0,0] | 0.654212
3 | [0,0,1,0,0] | 0.0947941
2 | [0,1,0,0,0] | 0.0946975
4 | [0,0,0,1,0] | 0.0904842
5 | [0,0,0,0,1] | 0.0221225
------------------------------------------------------------
------------------------------------------------------------
Component | Sobol index | Sobol total index
------------------------------------------------------------
0 | 0.662726 | 0.693362
1 | 0.0946975 | 0.10585
2 | 0.0947941 | 0.106069
3 | 0.0914871 | 0.10387
4 | 0.0221225 | 0.0253679
------------------------------------------------------------
draw Sobol' indices
.. code-block:: default
first_order = [sensitivityAnalysis.getSobolIndex(i) for i in range(dimension)]
total_order = [sensitivityAnalysis.getSobolTotalIndex(i) for i in range(dimension)]
graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(input_names, first_order, total_order)
view = viewer.View(graph)
.. image:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_functional_chaos_sensitivity_001.png
:alt: Sobol' indices
:class: sphx-glr-single-img
We saw that total order indices are close to first order,
so the higher order indices must be all quite close to 0
.. code-block:: default
for i in range(dimension):
for j in range(i):
print(input_names[i] + ' & '+ input_names[j], ":", sensitivityAnalysis.getSobolIndex([i, j]))
plt.show()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Hu & rw : 0.00939596043391626
Hl & rw : 0.0094957984784045
Hl & Hu : 0.0
L & rw : 0.00918479200468764
L & Hu : 0.0012912602896845825
L & Hl : 0.0013069732237138588
Kw & rw : 0.002220185764977052
Kw & Hu : 0.00031043066749530707
Kw & Hl : 0.0003119420529851785
Kw & L : 0.0003096614949037279
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 4.319 seconds)
.. _sphx_glr_download_auto_reliability_sensitivity_sensitivity_analysis_plot_functional_chaos_sensitivity.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_sensitivity.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_functional_chaos_sensitivity.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_