diff --git a/python/doc/examples/meta_modeling/kriging_metamodel/plot_kriging_cantilever_beam.py b/python/doc/examples/meta_modeling/kriging_metamodel/plot_kriging_cantilever_beam.py
index 5aee48dab43..d0c6b17d8df 100644
--- a/python/doc/examples/meta_modeling/kriging_metamodel/plot_kriging_cantilever_beam.py
+++ b/python/doc/examples/meta_modeling/kriging_metamodel/plot_kriging_cantilever_beam.py
@@ -132,7 +132,8 @@
 Y_test = model(X_test)
 
 # %%
-# The `MetaModelValidation` classe makes the validation easy. To create it, we use the validation samples and the metamodel.
+# The :class:`~openturns.MetaModelValidation` class makes surrogate model validation easy.
+# To create it, we use the validation samples and the metamodel.
 
 # %%
 val = ot.MetaModelValidation(X_test, Y_test, krigingMetamodel)
diff --git a/python/doc/examples/meta_modeling/kriging_metamodel/plot_kriging_cantilever_beam_hmat.py b/python/doc/examples/meta_modeling/kriging_metamodel/plot_kriging_cantilever_beam_hmat.py
index 005310bd00b..c63805cb111 100644
--- a/python/doc/examples/meta_modeling/kriging_metamodel/plot_kriging_cantilever_beam_hmat.py
+++ b/python/doc/examples/meta_modeling/kriging_metamodel/plot_kriging_cantilever_beam_hmat.py
@@ -144,7 +144,8 @@
 Y_test = model(X_test)
 
 # %%
-# The `MetaModelValidation` classe makes the validation easy. To create it, we use the validation samples and the metamodel.
+# The :class:`~openturns.MetaModelValidation` class is designed to validate the surrogate models.
+# To create it, we use a validation sample and a metamodel.
 
 # %%
 val = ot.MetaModelValidation(X_test, Y_test, krigingMetamodel)
diff --git a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_cv.py b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_cv.py
index 0af6c9977a2..0d26107d691 100644
--- a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_cv.py
+++ b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_cv.py
@@ -135,7 +135,70 @@ def compute_Q2_score_by_splitting(
 
 
 # %%
-# The next function computes the Q2 score by K-Fold.
+# The next function computes the mean squared error by K-Fold.
+
+# %%
+def computeMSENaiveKFold(
+    inputSample,
+    outputSample,
+    multivariateBasis,
+    totalDegree,
+    distribution,
+    kParameter=5,
+):
+    """
+    Compute mean squared error by (naive) KFold.
+
+    Parameters
+    ----------
+    inputSample : Sample(size, input_dimension)
+        The inputSample dataset.
+    outputSample : Sample(size, output_dimension)
+        The outputSample dataset.
+    multivariateBasis : multivariateBasis
+        The multivariate chaos multivariateBasis.
+    totalDegree : int
+        The total degree of the chaos polynomial.
+    distribution : Distribution.
+        The distribution of the input variable.
+    kParameter : int, in (2, sampleSize)
+        The parameter K.
+
+    Returns
+    -------
+    mse : Point(output_dimension)
+        The mean squared error.
+    """
+    #
+    sampleSize = inputSample.getSize()
+    outputDimension = outputSample.getDimension()
+    splitter = ot.KFoldSplitter(sampleSize, kParameter)
+    squaredResiduals = ot.Sample(sampleSize, outputDimension)
+    for indicesTrain, indicesTest in splitter:
+        inputSampleTrain, inputSampleTest = (
+            inputSample[indicesTrain],
+            inputSample[indicesTest],
+        )
+        outputSampleTrain, outputSampleTest = (
+            outputSample[indicesTrain],
+            outputSample[indicesTest],
+        )
+        chaosResultKFold = compute_sparse_least_squares_chaos(
+            inputSampleTrain,
+            outputSampleTrain,
+            multivariateBasis,
+            totalDegree,
+            distribution,
+        )
+        metamodelKFold = chaosResultKFold.getMetaModel()
+        predictionsKFold = metamodelKFold(inputSampleTest)
+        residualsKFold = outputSampleTest - predictionsKFold
+        foldSize = indicesTest.getSize()
+        for j in range(outputDimension):
+            for i in range(foldSize):
+                squaredResiduals[indicesTest[i], j] = residualsKFold[i, j] ** 2
+    mse = squaredResiduals.computeMean()
+    return mse
 
 
 # %%
diff --git a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_draw_validation.py b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_draw_validation.py
index 30831f02437..4eda73d5bd3 100644
--- a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_draw_validation.py
+++ b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_draw_validation.py
@@ -23,7 +23,8 @@
 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`.
+# The model 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
 input_names = im.distributionX.getDescription()