.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_meta_modeling/polynomial_chaos_metamodel/plot_chaos_sobol_confidence.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_meta_modeling_polynomial_chaos_metamodel_plot_chaos_sobol_confidence.py: Compute Sobol' indices confidence intervals ------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 7-32 In this example we will estimate confidence intervals for chaos Sobol' indices by bootstrap. Unlike estimation by sampling where the confidence intervals are obtained from the asymptotic distributions, there is currently no known analytical method to compute Sobol' confidence intervals, so we can fallback to bootstrap. Bootstraping the polynomial chaos is presented in [marelli2018]_ as "full bootstraping" and referred to as bootstrap-PCE or bPCE. Full bPCE can be CPU time consuming in some cases, e.g. when the dimension of the input random vector is large or when the training sample size is large. In the fast bPCE method, the sparse polynomial basis identified by the LARS algorithm is computed only once. Then bootstrapping is applied to the regression step on the sparse basis. Using fast bPCE is not straightforward, so we use full bPCE in the current example. To do this, we must simultaneously bootstrap in the input and output samples, so that the input/output mapping is preserved. In the example below, we show how to use the :class:`~openturns.BootstrapExperiment` class for this purpose. This involves the following steps: - Generate an initial input/output design of experiment. - Compute a fixed number of bootstrap samples from the original design. - For each input/output boostrap samples, compute Sobol' indices by functional chaos. - Compute quantiles of these Sobol' indices realizations to get the confidence intervals. .. GENERATED FROM PYTHON SOURCE LINES 34-38 .. code-block:: Python import openturns as ot from openturns.usecases import cantilever_beam import openturns.viewer as otv .. GENERATED FROM PYTHON SOURCE LINES 39-41 Load the cantilever beam use-case. We want to use the independent distribution to get meaningful Sobol' indices. .. GENERATED FROM PYTHON SOURCE LINES 41-46 .. code-block:: Python beam = cantilever_beam.CantileverBeam() g = beam.model distribution = beam.independentDistribution dim_input = distribution.getDimension() .. GENERATED FROM PYTHON SOURCE LINES 47-48 Create the tensorized polynomial basis, with the default linear enumeration function. .. GENERATED FROM PYTHON SOURCE LINES 48-51 .. code-block:: Python marginals = [distribution.getMarginal(i) for i in range(dim_input)] basis = ot.OrthogonalProductPolynomialFactory(marginals) .. GENERATED FROM PYTHON SOURCE LINES 52-53 Generate a learning sample with MC simulation (or retrieve the design from experimental data). .. GENERATED FROM PYTHON SOURCE LINES 53-58 .. code-block:: Python ot.RandomGenerator.SetSeed(0) N = 35 # size of the experimental design X = distribution.getSample(N) Y = g(X) .. GENERATED FROM PYTHON SOURCE LINES 59-106 .. code-block:: Python def computeSparseLeastSquaresChaos(X, Y, basis, total_degree, distribution): """ Create a sparse polynomial chaos with least squares. * Uses the enumeration rule from basis. * Uses LeastSquaresStrategy to compute the coefficients from linear least squares. * Uses LeastSquaresMetaModelSelectionFactory to select the polynomials in the basis using least angle regression stepwise (LARS) * Utilise LeastSquaresStrategy pour calculer les coefficients par moindres carrés. * Uses FixedStrategy to keep all coefficients that LARS has selected, up to the given maximum total degree. Parameters ---------- X : Sample(n) The input training design of experiments with n points Y : Sample(n) The input training design of experiments with n points basis : Basis The multivariate orthogonal polynomial basis total_degree : int The maximum total polynomial degree distribution : Distribution The distribution of the input random vector Returns ------- result : FunctionalChaosResult The polynomial chaos result """ selectionAlgorithm = ot.LeastSquaresMetaModelSelectionFactory() projectionStrategy = ot.LeastSquaresStrategy(selectionAlgorithm) enum_func = basis.getEnumerateFunction() P = enum_func.getBasisSizeFromTotalDegree(total_degree) adaptiveStrategy = ot.FixedStrategy(basis, P) algo = ot.FunctionalChaosAlgorithm( X, Y, distribution, adaptiveStrategy, projectionStrategy ) algo.run() result = algo.getResult() return result .. GENERATED FROM PYTHON SOURCE LINES 107-108 Build the chaos metamodel on the design of experiment. .. GENERATED FROM PYTHON SOURCE LINES 108-112 .. code-block:: Python total_degree = 3 result = computeSparseLeastSquaresChaos(X, Y, basis, total_degree, distribution) metamodel = result.getMetaModel() .. GENERATED FROM PYTHON SOURCE LINES 113-114 Generate a validation sample independent of the training sample. .. GENERATED FROM PYTHON SOURCE LINES 114-119 .. code-block:: Python n_valid = 1000 X_test = distribution.getSample(n_valid) Y_test = g(X_test) .. GENERATED FROM PYTHON SOURCE LINES 120-122 The MetaModelValidation class allows one to validate the metamodel on a test sample. Plot the observed versus the predicted outputs. .. GENERATED FROM PYTHON SOURCE LINES 122-129 .. code-block:: Python val = ot.MetaModelValidation(X_test, Y_test, metamodel) graph = val.drawValidation() Q2 = val.computePredictivityFactor()[0] graph.setTitle(f"Chaos validation - Q2={Q2*100.0:.2f}%") _ = otv.View(graph) .. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_sobol_confidence_001.png :alt: Chaos validation - Q2=99.90% :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_sobol_confidence_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 130-131 Retrieve Sobol' sensitivity measures associated to the polynomial chaos decomposition of the model. .. GENERATED FROM PYTHON SOURCE LINES 131-138 .. code-block:: Python chaosSI = ot.FunctionalChaosSobolIndices(result) fo_indices = [chaosSI.getSobolIndex(i) for i in range(dim_input)] to_indices = [chaosSI.getSobolTotalIndex(i) for i in range(dim_input)] input_names = g.getInputDescription() graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(input_names, fo_indices, to_indices) _ = otv.View(graph) .. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_sobol_confidence_002.png :alt: Sobol' indices :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_sobol_confidence_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 139-142 We define the `multiBootstrap` function in order to simultaneously bootstrap in the input and output samples, so that the input/output mapping is preserved. We use the `GenerateSelection` method of the :class:`~openturns.BootstrapExperiment` class. .. GENERATED FROM PYTHON SOURCE LINES 142-164 .. code-block:: Python def multiBootstrap(*data): """ Bootstrap multiple samples at once. Parameters ---------- data : sequence of Sample Multiple samples to bootstrap. Returns ------- data_boot : sequence of Sample The bootstrap samples. """ assert len(data) > 0, "empty list" size = data[0].getSize() selection = ot.BootstrapExperiment.GenerateSelection(size, size) return [Z[selection] for Z in data] .. GENERATED FROM PYTHON SOURCE LINES 165-166 Generate input/output bootstrap samples from the initial design. .. GENERATED FROM PYTHON SOURCE LINES 166-170 .. code-block:: Python X_boot, Y_boot = multiBootstrap(X, Y) print(X_boot[:5]) print(Y_boot[:5]) .. rst-class:: sphx-glr-script-out .. code-block:: none [ E F L I ] 0 : [ 7.09376e+10 305.568 2.55637 1.34556e-07 ] 1 : [ 6.5963e+10 228.332 2.57804 1.53991e-07 ] 2 : [ 6.60008e+10 325.25 2.5146 1.45825e-07 ] 3 : [ 6.52443e+10 320.71 2.54539 1.40993e-07 ] 4 : [ 6.86026e+10 282.187 2.59881 1.42475e-07 ] [ y0 ] 0 : [ 0.178271 ] 1 : [ 0.128386 ] 2 : [ 0.179111 ] 3 : [ 0.191651 ] 4 : [ 0.168912 ] .. GENERATED FROM PYTHON SOURCE LINES 171-201 .. code-block:: Python def computeChaosSensitivity(X, Y, basis, total_degree, distribution): """ Compute the first and total order Sobol' indices from a polynomial chaos. Parameters ---------- X, Y : Sample Input/Output design basis : Basis Tensorized basis total_degree : int Maximal total degree distribution : Distribution Input distribution Returns ------- first_order, total_order: list of float The first and total order indices. """ dim_input = X.getDimension() result = computeSparseLeastSquaresChaos(X, Y, basis, total_degree, distribution) chaosSI = ot.FunctionalChaosSobolIndices(result) first_order = [chaosSI.getSobolIndex(i) for i in range(dim_input)] total_order = [chaosSI.getSobolTotalIndex(i) for i in range(dim_input)] return first_order, total_order .. GENERATED FROM PYTHON SOURCE LINES 202-238 .. code-block:: Python def computeBootstrapChaosSobolIndices( X, Y, basis, total_degree, distribution, bootstrap_size ): """ Computes a bootstrap sample of first and total order indices from polynomial chaos. Parameters ---------- X, Y : Sample Input/Output design basis : Basis Tensorized basis total_degree : int Maximal total degree distribution : Distribution Input distribution bootstrap_size : interval The bootstrap sample size """ dim_input = X.getDimension() fo_sample = ot.Sample(0, dim_input) to_sample = ot.Sample(0, dim_input) unit_eps = ot.Interval([1e-9] * dim_input, [1 - 1e-9] * dim_input) for i in range(bootstrap_size): X_boot, Y_boot = multiBootstrap(X, Y) first_order, total_order = computeChaosSensitivity( X_boot, Y_boot, basis, total_degree, distribution ) if unit_eps.contains(first_order) and unit_eps.contains(total_order): fo_sample.add(first_order) to_sample.add(total_order) return fo_sample, to_sample .. GENERATED FROM PYTHON SOURCE LINES 239-240 Compute sample of Sobol' indices by boostrap. .. GENERATED FROM PYTHON SOURCE LINES 240-245 .. code-block:: Python bootstrap_size = 500 fo_sample, to_sample = computeBootstrapChaosSobolIndices( X, Y, basis, total_degree, distribution, bootstrap_size ) .. GENERATED FROM PYTHON SOURCE LINES 246-293 .. code-block:: Python def computeSobolIndicesConfidenceInterval(fo_sample, to_sample, alpha=0.95): """ From a sample of first or total order indices, compute a bilateral confidence interval of level alpha. Estimates the distribution of the first and total order Sobol' indices from a bootstrap estimation. Then computes a bilateral confidence interval for each marginal. Parameters ---------- fo_sample: ot.Sample(n, dim_input) The first order indices to_sample: ot.Sample(n, dim_input) The total order indices alpha : float The confidence level Returns ------- fo_interval : ot.Interval The confidence interval of first order Sobol' indices to_interval : ot.Interval The confidence interval of total order Sobol' indices """ dim_input = fo_sample.getDimension() fo_lb = [0] * dim_input fo_ub = [0] * dim_input to_lb = [0] * dim_input to_ub = [0] * dim_input for i in range(dim_input): fo_i = fo_sample[:, i] to_i = to_sample[:, i] beta = (1.0 - alpha) / 2 fo_lb[i] = fo_i.computeQuantile(beta)[0] fo_ub[i] = fo_i.computeQuantile(1.0 - beta)[0] to_lb[i] = to_i.computeQuantile(beta)[0] to_ub[i] = to_i.computeQuantile(1.0 - beta)[0] # Create intervals fo_interval = ot.Interval(fo_lb, fo_ub) to_interval = ot.Interval(to_lb, to_ub) return fo_interval, to_interval .. GENERATED FROM PYTHON SOURCE LINES 294-295 Compute confidence intervals from the indices samples. .. GENERATED FROM PYTHON SOURCE LINES 295-297 .. code-block:: Python fo_interval, to_interval = computeSobolIndicesConfidenceInterval(fo_sample, to_sample) .. GENERATED FROM PYTHON SOURCE LINES 298-330 .. code-block:: Python def computeAndDrawSobolIndices( N, basis, total_degree, distribution, bootstrap_size=500, alpha=0.95 ): """ Compute and draw Sobol' indices from a polynomial chaos based on a given sample size. Compute confidence intervals at level alpha from bootstrap. """ X = distribution.getSample(N) Y = g(X) fo_sample, to_sample = computeBootstrapChaosSobolIndices( X, Y, basis, total_degree, distribution, bootstrap_size ) fo_interval, to_interval = computeSobolIndicesConfidenceInterval( fo_sample, to_sample, alpha ) graph = ot.SobolIndicesAlgorithm.DrawSobolIndices( input_names, fo_sample.computeMean(), to_sample.computeMean(), fo_interval, to_interval, ) graph.setTitle(f"Sobol' indices - N={N}") graph.setIntegerXTick(True) return graph .. GENERATED FROM PYTHON SOURCE LINES 331-333 Plot the boostrap Sobol' indices mean and the confidence intervals. The confidence length will shrink if we increase the initial sample size. .. GENERATED FROM PYTHON SOURCE LINES 335-336 sphinx_gallery_thumbnail_number = 3 .. GENERATED FROM PYTHON SOURCE LINES 336-341 .. code-block:: Python graph = computeAndDrawSobolIndices(30, basis, total_degree, distribution) _ = otv.View(graph) otv.View.ShowAll() .. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_sobol_confidence_003.png :alt: Sobol' indices - N=30 :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_sobol_confidence_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 342-349 We see that the variable F has the highest sensitivity indices and that confidence intervals do not change this conclusion. The confidence intervals of the sensitivity indices of the variables L and I are similar so that we cannot say which of L or I is more significant than the other: both variables have similar sensitivity indices. The least sensitive variable is E, but the confidence intervals do not cross the X axis. Hence, there is no evidence that the Sobol' indices of E are zero. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.766 seconds) .. _sphx_glr_download_auto_meta_modeling_polynomial_chaos_metamodel_plot_chaos_sobol_confidence.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_chaos_sobol_confidence.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_chaos_sobol_confidence.py `