.. 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_cv.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_cv.py: Polynomial chaos expansion cross-validation =========================================== .. GENERATED FROM PYTHON SOURCE LINES 7-34 Introduction ------------ In this example, we show how to perform the cross-validation of the :ref:`Ishigami model` using polynomial chaos expansion. More precisely, we use the methods suggested in [muller2016]_, chapter 5, page 257. We make the assumption that a dataset is given and we create a metamodel using this data. Once done, we want to assess the quality of the metamodel using the data we have. Another example of this method is presented in :doc:`/auto_numerical_methods/general_methods/plot_pce_design`. In this example, we compare two methods: - split the data into two subsets, one for training and one for testing, - use k-fold validation. The split of the data is performed by the `compute_Q2_score_by_splitting` function below. It uses 75% of the data to estimate the coefficients of the metamodel (this is the training step) and use 25% of the data to estimate the :math:`Q^2` score (this is the validation step). To do this, we use the `split` method of the :class:`~openturns.Sample`. The K-Fold validation is performed by the `compute_Q2_score_by_kfold` function below. It uses the K-Fold method with :math:`k = 5`. The code uses the :class:`~openturns.KFoldSplitter` class, which computes the appropriate indices. Similar results can be obtained with :class:`~openturns.LeaveOneOutSplitter` at a higher cost. This cross-validation method is used here to see which polynomial degree leads to an accurate metamodel of the Ishigami test function. .. GENERATED FROM PYTHON SOURCE LINES 36-40 .. code-block:: Python import openturns as ot import openturns.viewer as otv from openturns.usecases import ishigami_function .. GENERATED FROM PYTHON SOURCE LINES 41-43 Define helper functions ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 45-47 The next function creates a polynomial chaos expansion from a training data set and a total degree. .. GENERATED FROM PYTHON SOURCE LINES 49-97 .. code-block:: Python def compute_sparse_least_squares_chaos( inputTrain, outputTrain, basis, totalDegree, distribution ): """ Create a sparse polynomial chaos based on least squares. * Uses the enumerate rule in basis. * Uses the LeastSquaresStrategy to compute the coefficients based on least squares. * Uses LeastSquaresMetaModelSelectionFactory to use the LARS selection method. * Uses FixedStrategy in order to keep all the coefficients that the LARS method selected. Parameters ---------- inputTrain : Sample The input design of experiments. outputTrain : Sample The output design of experiments. basis : Basis The multivariate chaos basis. totalDegree : int The total degree of the chaos polynomial. distribution : Distribution. The distribution of the input variable. Returns ------- result : PolynomialChaosResult The estimated polynomial chaos. """ selectionAlgorithm = ot.LeastSquaresMetaModelSelectionFactory() projectionStrategy = ot.LeastSquaresStrategy( inputTrain, outputTrain, selectionAlgorithm ) enumerateFunction = basis.getEnumerateFunction() basisSize = enumerateFunction.getBasisSizeFromTotalDegree(totalDegree) adaptiveStrategy = ot.FixedStrategy(basis, basisSize) chaosAlgo = ot.FunctionalChaosAlgorithm( inputTrain, outputTrain, distribution, adaptiveStrategy, projectionStrategy ) chaosAlgo.run() result = chaosAlgo.getResult() return result .. GENERATED FROM PYTHON SOURCE LINES 98-100 The next function computes the Q2 score by splitting the data set into a training set and a test set. .. GENERATED FROM PYTHON SOURCE LINES 103-137 .. code-block:: Python def compute_Q2_score_by_splitting( X, Y, basis, totalDegree, distribution, split_fraction=0.75 ): """ Compute Q2 score by splitting into train/test sets. Parameters ---------- X : Sample(size, input_dimension) The X dataset. Y : Sample(size, output_dimension) The Y dataset. Returns ------- Q2_score : float The Q2 score. """ training_sample_size = X.getSize() X_train = ot.Sample(X) Y_train = ot.Sample(Y) split_index = int(split_fraction * training_sample_size) X_test = X_train.split(split_index) Y_test = Y_train.split(split_index) result = compute_sparse_least_squares_chaos( X_train, Y_train, basis, totalDegree, distribution ) metamodel = result.getMetaModel() val = ot.MetaModelValidation(X_test, Y_test, metamodel) Q2_score = val.computePredictivityFactor()[0] return Q2_score .. GENERATED FROM PYTHON SOURCE LINES 138-139 The next function computes the Q2 score by K-Fold. .. GENERATED FROM PYTHON SOURCE LINES 142-176 .. code-block:: Python def compute_Q2_score_by_kfold(X, Y, basis, totalDegree, distribution, n_folds=5): """ Compute score by KFold. Parameters ---------- X : Sample(size, input_dimension) The X dataset. Y : Sample(size, output_dimension) The Y dataset. Returns ------- Q2_score : float The Q2 score. """ # training_sample_size = X.getSize() splitter = ot.KFoldSplitter(training_sample_size, n_folds) Q2_score_list = ot.Sample(0, 1) for indices1, indices2 in splitter: X_train, X_test = X[indices1], X[indices2] Y_train, Y_test = Y[indices1], Y[indices2] result = compute_sparse_least_squares_chaos( X_train, Y_train, basis, totalDegree, distribution ) metamodel = result.getMetaModel() val = ot.MetaModelValidation(X_test, Y_test, metamodel) Q2_local = val.computePredictivityFactor()[0] Q2_score_list.add([Q2_local]) Q2_score = Q2_score_list.computeMean()[0] return Q2_score .. GENERATED FROM PYTHON SOURCE LINES 177-179 Define the training data set ---------------------------- .. GENERATED FROM PYTHON SOURCE LINES 181-182 We start by generating the input and output sample. We use a sample size equal to 1000. .. GENERATED FROM PYTHON SOURCE LINES 184-193 .. code-block:: Python im = ishigami_function.IshigamiModel() im.distributionX.setDescription(["X0", "X1", "X2"]) im.model.setOutputDescription(["$Y$"]) ot.RandomGenerator.SetSeed(0) sample_size = 500 X = im.distributionX.getSample(sample_size) print("Input sample:") X[:5] .. rst-class:: sphx-glr-script-out .. code-block:: none Input sample: .. raw:: html
X0X1X2
00.8160385-1.5469512.67553
12.405236-0.4441351-2.192122
2-2.291626-3.0132041.30032
3-2.9373720.0734501-0.7336488
4-0.960969-0.1161257-0.4336099


.. GENERATED FROM PYTHON SOURCE LINES 194-200 .. code-block:: Python Y = im.model(X) Y.setDescription(["Y"]) print("Output sample:") Y[:5] .. rst-class:: sphx-glr-script-out .. code-block:: none Output sample: .. raw:: html
Y
011.45723
13.514782
2-0.8512832
3-0.170983
4-0.728672


.. GENERATED FROM PYTHON SOURCE LINES 201-204 We compute a polynomial chaos decomposition with a total degree equal to 5. In order to reduce the number of coefficients, the `compute_sparse_least_squares_chaos` function creates a sparse chaos using regression and the LARS method. .. GENERATED FROM PYTHON SOURCE LINES 206-214 .. code-block:: Python dimension = im.distributionX.getDimension() basis = ot.OrthogonalProductPolynomialFactory( [im.distributionX.getMarginal(i) for i in range(dimension)] ) totalDegree = 5 # Polynomial degree result = compute_sparse_least_squares_chaos(X, Y, basis, totalDegree, im.distributionX) result .. raw:: html
FunctionalChaosResult
  • input dimension: 3
  • output dimension: 1
  • distribution dimension: 3
  • transformation: 3 -> 3
  • inverse transformation: 3 -> 3
  • orthogonal basis dimension: 3
  • indices size: 12
  • relative errors: [0.000326285]
  • residuals: [0.0644485]
Index Multi-index Coeff.
0 [0,0,0] 3.661918
1 [1,0,0] 1.599859
2 [1,0,1] -0.1517826
3 [0,2,0] -0.616015
4 [0,0,2] -0.1429231
5 [3,0,0] -1.21482
6 [1,0,2] 1.421431
7 [0,4,0] -1.983213
8 [5,0,0] 0.161996
9 [4,1,0] 0.1837313
10 [3,0,2] -1.165708
11 [1,0,4] 0.3744418


.. GENERATED FROM PYTHON SOURCE LINES 215-216 Get the metamodel. .. GENERATED FROM PYTHON SOURCE LINES 218-220 .. code-block:: Python metamodel = result.getMetaModel() .. GENERATED FROM PYTHON SOURCE LINES 221-223 Validate the metamodel from a test set -------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 225-227 In order to validate our polynomial chaos, we generate an extra pair of input / output samples and use the :class:`~openturns.MetaModelValidation` class. .. GENERATED FROM PYTHON SOURCE LINES 227-237 .. code-block:: Python test_sample_size = 200 # Size of the validation design of experiments inputTest = im.distributionX.getSample(test_sample_size) outputTest = im.model(inputTest) validation = ot.MetaModelValidation(inputTest, outputTest, metamodel) Q2 = validation.computePredictivityFactor()[0] graph = validation.drawValidation() graph.setTitle("Q2=%.2f, n=%d" % (Q2, test_sample_size)) view = otv.View(graph) .. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_cv_001.png :alt: Q2=0.84, n=200 :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_cv_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 238-256 The plot shows that the score is relatively high and might be satisfactory for some purposes. There are however several issues with the previous procedure: - It may happen that the data in the validation sample with size 200 is more difficult to fit than the data in the training dataset. In this case, the :math:`Q^2` score may be pessimistic. - It may happen that the data in the validation sample with size 200 is less difficult to fit than the data in the validation dataset. In this case, the :math:`Q^2` score may be optimistic. - We may not be able to generate an extra dataset for validation. In this case, a part of the original dataset should be used for validation. - The polynomial degree may not be appropriate for this data. - The dataset may be ordered in some way, so that the split is somewhat arbitrary. One solution to circumvent this is to randomly shuffle the data. This can be done using the :class:`~openturns.KPermutationsDistribution`. The K-Fold validation aims at solving some of these issues, so that all the available data is used in order to estimate the :math:`Q^2` score. .. GENERATED FROM PYTHON SOURCE LINES 258-260 Compute the Q2 score from a test set ------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 262-263 In the following script, we compute the :math:`Q^2` score associated with each polynomial degree from 1 to 10. .. GENERATED FROM PYTHON SOURCE LINES 263-275 .. code-block:: Python degree_max = 10 degree_list = list(range(1, degree_max)) n_degrees = len(degree_list) score_sample = ot.Sample(len(degree_list), 1) for i in range(n_degrees): totalDegree = degree_list[i] score_sample[i, 0] = compute_Q2_score_by_splitting( X, Y, basis, totalDegree, im.distributionX ) print(f"split - degree = {totalDegree}, score = {score_sample[i, 0]:.4f}") .. rst-class:: sphx-glr-script-out .. code-block:: none split - degree = 1, score = 0.1861 split - degree = 2, score = 0.1861 split - degree = 3, score = 0.4948 split - degree = 4, score = 0.7338 split - degree = 5, score = 0.8272 split - degree = 6, score = 0.9772 split - degree = 7, score = 0.9889 split - degree = 8, score = 0.9994 split - degree = 9, score = 0.9997 .. GENERATED FROM PYTHON SOURCE LINES 276-281 .. code-block:: Python graph = ot.Graph("Split", "Degree", "$Q^2$", True) cloud = ot.Cloud(ot.Sample.BuildFromPoint(degree_list), score_sample) graph.add(cloud) view = otv.View(graph) .. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_cv_002.png :alt: Split :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_cv_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 282-284 We see that the polynomial degree may be increased up to degree 7, after which the :math:`Q^2` score does not increase much. .. GENERATED FROM PYTHON SOURCE LINES 286-291 Compute the Q2 score from K-Fold cross-validation ------------------------------------------------- One limitation of the previous method is that the estimate of the :math:`Q^2` may be sensitive to the particular split of the dataset. The following script uses 5-Fold cross validation to estimate the :math:`Q^2` score. .. GENERATED FROM PYTHON SOURCE LINES 293-301 .. code-block:: Python score_sample = ot.Sample(len(degree_list), 1) for i in range(n_degrees): totalDegree = degree_list[i] score_sample[i, 0] = compute_Q2_score_by_kfold( X, Y, basis, totalDegree, im.distributionX ) print(f"k-fold, degree = {totalDegree}, score = {score_sample[i, 0]:.4f}") .. rst-class:: sphx-glr-script-out .. code-block:: none k-fold, degree = 1, score = 0.1516 k-fold, degree = 2, score = 0.1499 k-fold, degree = 3, score = 0.4480 k-fold, degree = 4, score = 0.7120 k-fold, degree = 5, score = 0.8245 k-fold, degree = 6, score = 0.9805 k-fold, degree = 7, score = 0.9901 k-fold, degree = 8, score = 0.9996 k-fold, degree = 9, score = 0.9998 .. GENERATED FROM PYTHON SOURCE LINES 302-307 .. code-block:: Python graph = ot.Graph("K-Fold", "Degree", "$Q^2$", True) cloud = ot.Cloud(ot.Sample.BuildFromPoint(degree_list), score_sample) graph.add(cloud) view = otv.View(graph) .. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_cv_003.png :alt: K-Fold :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_chaos_cv_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 308-309 The conclusion is similar to the previous method. .. GENERATED FROM PYTHON SOURCE LINES 311-325 Conclusion ---------- When we select the best polynomial degree which maximizes the :math:`Q^2` score, the danger is that the validation set is used both for computing the :math:`Q^2` and to maximize it: hence, the :math:`Q^2` score may be optimistic. In [muller2016]_, chapter 5, page 269, the authors advocate the split of the dataset into three subsets: - the training set, - the validation set, - the test set. When selecting the best parameters, the validation set is used. When estimating the :math:`Q^2` score with the best parameters, the test set is used. .. GENERATED FROM PYTHON SOURCE LINES 327-328 .. code-block:: Python otv.View.ShowAll() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 11.031 seconds) .. _sphx_glr_download_auto_meta_modeling_polynomial_chaos_metamodel_plot_chaos_cv.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_cv.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_chaos_cv.py `