.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_meta_modeling/kriging_metamodel/plot_kriging_categorical.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_kriging_metamodel_plot_kriging_categorical.py: Kriging: metamodel with continuous and categorical variables ============================================================ .. GENERATED FROM PYTHON SOURCE LINES 7-10 We consider here the surrogate modeling of an analytical function characterized by continuous and categorical variables .. GENERATED FROM PYTHON SOURCE LINES 12-20 .. code-block:: Python import openturns as ot import openturns.experimental as otexp import numpy as np import matplotlib.pyplot as plt # Seed chosen in order to obtain a visually nice plot ot.RandomGenerator.SetSeed(5) .. GENERATED FROM PYTHON SOURCE LINES 21-24 We first show the advantage of modeling the various levels of a mixed continuous / categorical function through a single surrogate model on a simple test-case taken from [pelamatti2020]_, defined below. .. GENERATED FROM PYTHON SOURCE LINES 27-41 .. code-block:: Python def illustrativeFunc(inp): x, z = inp y = np.cos(7 * x) + 0.5 * z return [y] dim = 2 fun = ot.PythonFunction(dim, 1, illustrativeFunc) numberOfZLevels = 2 # Number of categorical levels for z # Input distribution dist = ot.JointDistribution( [ot.Uniform(0, 1), ot.UserDefined(ot.Sample.BuildFromPoint(range(numberOfZLevels)))] ) .. GENERATED FROM PYTHON SOURCE LINES 42-45 In this example, we compare the performances of the :class:`~openturns.experimental.LatentVariableModel` with a naive approach, which would consist in modeling each combination of categorical variables through a separate and independent Gaussian process. .. GENERATED FROM PYTHON SOURCE LINES 47-51 In order to deal with mixed continuous / categorical problems we can rely on the :class:`~openturns.ProductCovarianceModel` class. We start here by defining the product kernel, which combines :class:`~openturns.SquaredExponential` kernels for the continuous variables, and :class:`~openturns.experimental.LatentVariableModel` for the categorical ones. .. GENERATED FROM PYTHON SOURCE LINES 53-71 .. code-block:: Python latDim = 1 # Dimension of the latent space activeCoord = 1 + latDim * ( numberOfZLevels - 2 ) # Nb of active coordinates in the latent space kx = ot.SquaredExponential(1) kz = otexp.LatentVariableModel(numberOfZLevels, latDim) kLV = ot.ProductCovarianceModel([kx, kz]) kLV.setNuggetFactor(1e-6) # Bounds for the hyperparameter optimization lowerBoundLV = [1e-4] * dim + [-10.0] * activeCoord upperBoundLV = [2.0] * dim + [10.0] * activeCoord boundsLV = ot.Interval(lowerBoundLV, upperBoundLV) # Distribution for the hyperparameters initialization initDistLV = ot.DistributionCollection() for i in range(len(lowerBoundLV)): initDistLV.add(ot.Uniform(lowerBoundLV[i], upperBoundLV[i])) initDistLV = ot.JointDistribution(initDistLV) .. GENERATED FROM PYTHON SOURCE LINES 72-74 As a reference, we consider a purely continuous kernel for independent Gaussian processes. One for each combination of categorical variables levels. .. GENERATED FROM PYTHON SOURCE LINES 76-87 .. code-block:: Python kIndependent = ot.SquaredExponential(1) lowerBoundInd = [1e-4] upperBoundInd = [20.0] boundsInd = ot.Interval(lowerBoundInd, upperBoundInd) initDistInd = ot.DistributionCollection() for i in range(len(lowerBoundInd)): initDistInd.add(ot.Uniform(lowerBoundInd[i], upperBoundInd[i])) initDistInd = ot.JointDistribution(initDistInd) initSampleInd = initDistInd.getSample(10) optAlgInd = ot.MultiStart(ot.Cobyla(), initSampleInd) .. GENERATED FROM PYTHON SOURCE LINES 88-89 Generate the training data set .. GENERATED FROM PYTHON SOURCE LINES 89-97 .. code-block:: Python x = dist.getSample(10) y = fun(x) # And the plotting data set xPlt = dist.getSample(200) xPlt = xPlt.sort() yPlt = fun(xPlt) .. GENERATED FROM PYTHON SOURCE LINES 98-99 Initialize and parameterize the optimization algorithm .. GENERATED FROM PYTHON SOURCE LINES 99-102 .. code-block:: Python initSampleLV = initDistLV.getSample(30) optAlgLV = ot.MultiStart(ot.Cobyla(), initSampleLV) .. GENERATED FROM PYTHON SOURCE LINES 103-104 Create and train the Gaussian process models .. GENERATED FROM PYTHON SOURCE LINES 104-127 .. code-block:: Python basis = ot.ConstantBasisFactory(2).build() algoLV = ot.KrigingAlgorithm(x, y, kLV, basis) algoLV.setOptimizationAlgorithm(optAlgLV) algoLV.setOptimizationBounds(boundsLV) algoLV.run() resLV = algoLV.getResult() algoIndependentList = [] for z in range(2): # Select the training samples corresponding to the correct combination # of categorical levels ind = np.where(np.all(np.array(x[:, 1]) == z, axis=1))[0] xLoc = x[ind][:, 0] yLoc = y[ind] # Create and train the Gaussian process models basis = ot.ConstantBasisFactory(1).build() algoIndependent = ot.KrigingAlgorithm(xLoc, yLoc, kIndependent, basis) algoIndependent.setOptimizationAlgorithm(optAlgInd) algoIndependent.setOptimizationBounds(boundsInd) algoIndependent.run() algoIndependentList.append(algoIndependent.getResult()) .. GENERATED FROM PYTHON SOURCE LINES 128-130 Plot the prediction of the mixed continuous / categorical GP, as well as the one of the two separate continuous GPs .. GENERATED FROM PYTHON SOURCE LINES 130-183 .. code-block:: Python fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(15, 10)) for z in range(numberOfZLevels): # Select the training samples corresponding to the correct combination # of categorical levels ind = np.where(np.all(np.array(x[:, 1]) == z, axis=1))[0] xLoc = x[ind][:, 0] yLoc = y[ind] # Compute the models predictive performances on a validation data set. # The predictions are computed independently for each level of z, # i.e., by only considering the values of z corresponding to the # target level. ind = np.where(np.all(np.array(xPlt[:, 1]) == z, axis=1))[0] xPltInd = xPlt[ind] yPltInd = yPlt[ind] predMeanLV = resLV.getConditionalMean(xPltInd) predMeanInd = algoIndependentList[z].getConditionalMean(xPltInd[:, 0]) predSTDLV = np.sqrt(resLV.getConditionalMarginalVariance(xPltInd)) predSTDInd = np.sqrt( algoIndependentList[z].getConditionalMarginalVariance(xPltInd[:, 0]) ) (trainingData,) = ax1.plot(xLoc[:, 0], yLoc, "r*") (trueFunction,) = ax1.plot(xPltInd[:, 0], yPltInd, "k--") (prediction,) = ax1.plot(xPltInd[:, 0], predMeanLV, "b-") stdPred = ax1.fill_between( xPltInd[:, 0].asPoint(), (predMeanLV - predSTDLV).asPoint(), (predMeanLV + predSTDLV).asPoint(), alpha=0.5, color="blue", ) ax2.plot(xLoc[:, 0], yLoc, "r*") ax2.plot(xPltInd[:, 0], yPltInd, "k--") ax2.plot(xPltInd[:, 0], predMeanInd, "b-") ax2.fill_between( xPltInd[:, 0].asPoint(), (predMeanInd - predSTDInd).asPoint(), (predMeanInd + predSTDInd).asPoint(), alpha=0.5, color="blue", ) ax1.legend( [trainingData, trueFunction, prediction, stdPred], ["Training data", "True function", "Prediction", "Prediction standard deviation"], ) ax1.set_title("Mixed continuous-categorical modeling") ax2.set_title("Separate modeling") ax2.set_xlabel("x", fontsize=15) ax1.set_ylabel("y", fontsize=15) ax2.set_ylabel("y", fontsize=15) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_categorical_001.png :alt: Mixed continuous-categorical modeling, Separate modeling :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_categorical_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(138.47222222222223, 0.5, 'y') .. GENERATED FROM PYTHON SOURCE LINES 184-187 It can be seen that the joint modeling of categorical and continuous variables improves the overall prediction accuracy, as the Gaussian process model is able to exploit the information provided by the entire training data set. .. GENERATED FROM PYTHON SOURCE LINES 189-192 We now consider a more complex function which is a modified version of the Goldstein function, taken from [pelamatti2020]_. This function depends on 2 continuous variables and 2 categorical ones. Each categorical variable is characterized by 3 levels. .. GENERATED FROM PYTHON SOURCE LINES 195-259 .. code-block:: Python def h(x1, x2, x3, x4): y = ( 53.3108 + 0.184901 * x1 - 5.02914 * x1**3 * 1e-6 + 7.72522 * x1**4 * 1e-8 - 0.0870775 * x2 - 0.106959 * x3 + 7.98772 * x3**3 * 1e-6 + 0.00242482 * x4 + 1.32851 * x4**3 * 1e-6 * 0.00146393 * x1 * x2 - 0.00301588 * x1 * x3 - 0.00272291 * x1 * x4 + 0.0017004 * x2 * x3 + 0.0038428 * x2 * x4 - 0.000198969 * x3 * x4 + 1.86025 * x1 * x2 * x3 * 1e-5 - 1.88719 * x1 * x2 * x4 * 1e-6 + 2.50923 * x1 * x3 * x4 * 1e-5 - 5.62199 * x2 * x3 * x4 * 1e-5 ) return y def Goldstein(inp): x1, x2, z1, z2 = inp x1 = 100 * x1 x2 = 100 * x2 if z1 == 0: x3 = 80 elif z1 == 1: x3 = 20 elif z1 == 2: x3 = 50 else: print("error, no matching category z1") if z2 == 0: x4 = 20 elif z2 == 1: x4 = 80 elif z2 == 2: x4 = 50 else: print("error, no matching category z2") return [h(x1, x2, x3, x4)] dim = 4 fun = ot.PythonFunction(dim, 1, Goldstein) numberOfZLevels1 = 3 # Number of categorical levels for z1 numberOfZLevels2 = 3 # Number of categorical levels for z2 # Input distribution dist = ot.JointDistribution( [ ot.Uniform(0, 1), ot.Uniform(0, 1), ot.UserDefined(ot.Sample.BuildFromPoint(range(numberOfZLevels1))), ot.UserDefined(ot.Sample.BuildFromPoint(range(numberOfZLevels2))), ] ) .. GENERATED FROM PYTHON SOURCE LINES 260-263 As in the previous example, we start here by defining the product kernel, which combines :class:`~openturns.SquaredExponential` kernels for the continuous variables, and :class:`~openturns.experimental.LatentVariableModel` for the categorical ones. .. GENERATED FROM PYTHON SOURCE LINES 265-285 .. code-block:: Python latDim = 2 # Dimension of the latent space activeCoord = ( 2 + latDim * (numberOfZLevels1 - 2) + latDim * (numberOfZLevels2 - 2) ) # Nb ative coordinates in the latent space kx1 = ot.SquaredExponential(1) kx2 = ot.SquaredExponential(1) kz1 = otexp.LatentVariableModel(numberOfZLevels1, latDim) kz2 = otexp.LatentVariableModel(numberOfZLevels2, latDim) kLV = ot.ProductCovarianceModel([kx1, kx2, kz1, kz2]) kLV.setNuggetFactor(1e-6) # Bounds for the hyperparameter optimization lowerBoundLV = [1e-4] * dim + [-10] * activeCoord upperBoundLV = [3.0] * dim + [10.0] * activeCoord boundsLV = ot.Interval(lowerBoundLV, upperBoundLV) # Distribution for the hyperparameters initialization initDistLV = ot.DistributionCollection() for i in range(len(lowerBoundLV)): initDistLV.add(ot.Uniform(lowerBoundLV[i], upperBoundLV[i])) initDistLV = ot.JointDistribution(initDistLV) .. GENERATED FROM PYTHON SOURCE LINES 286-288 Alternatively, we consider a purely continuous kernel for independent Gaussian processes. one for each combination of categorical variables levels. .. GENERATED FROM PYTHON SOURCE LINES 290-301 .. code-block:: Python kIndependent = ot.SquaredExponential(2) lowerBoundInd = [1e-4, 1e-4] upperBoundInd = [3.0, 3.0] boundsInd = ot.Interval(lowerBoundInd, upperBoundInd) initDistInd = ot.DistributionCollection() for i in range(len(lowerBoundInd)): initDistInd.add(ot.Uniform(lowerBoundInd[i], upperBoundInd[i])) initDistInd = ot.JointDistribution(initDistInd) initSampleInd = initDistInd.getSample(10) optAlgInd = ot.MultiStart(ot.Cobyla(), initSampleInd) .. GENERATED FROM PYTHON SOURCE LINES 302-306 In order to assess their respective robustness with regards to the training data set, we repeat the experiments 3 times with different training of size 72, and compute each time the normalized prediction Root Mean Squared Error (RMSE) on a test data set of size 1000. .. GENERATED FROM PYTHON SOURCE LINES 306-368 .. code-block:: Python rmseLVList = [] rmseIndList = [] for rep in range(3): # Generate the normalized training data set x = dist.getSample(72) y = fun(x) yMax = y.getMax() yMin = y.getMin() y = (y - yMin) / (yMin - yMax) # Initialize and parameterize the optimization algorithm initSampleLV = initDistLV.getSample(10) optAlgLV = ot.MultiStart(ot.Cobyla(), initSampleLV) # Create and train the Gaussian process models basis = ot.ConstantBasisFactory(dim).build() algoLV = ot.KrigingAlgorithm(x, y, kLV, basis) algoLV.setOptimizationAlgorithm(optAlgLV) algoLV.setOptimizationBounds(boundsLV) algoLV.run() resLV = algoLV.getResult() # Compute the models predictive performances on a validation data set xVal = dist.getSample(1000) yVal = fun(xVal) yVal = (yVal - yMin) / (yMin - yMax) valLV = ot.MetaModelValidation(yVal, resLV.getMetaModel()(xVal)) rmseLV = valLV.getResidualSample().computeStandardDeviation()[0] rmseLVList.append(rmseLV) error = ot.Sample(0, 1) for z1 in range(numberOfZLevels1): for z2 in range(numberOfZLevels2): # Select the training samples corresponding to the correct combination # of categorical levels ind = np.where(np.all(np.array(x[:, 2:]) == [z1, z2], axis=1))[0] xLoc = x[ind][:, :2] yLoc = y[ind] # Create and train the Gaussian process models basis = ot.ConstantBasisFactory(2).build() algoIndependent = ot.KrigingAlgorithm(xLoc, yLoc, kIndependent, basis) algoIndependent.setOptimizationAlgorithm(optAlgInd) algoIndependent.setOptimizationBounds(boundsInd) algoIndependent.run() resInd = algoIndependent.getResult() # Compute the models predictive performances on a validation data set ind = np.where(np.all(np.array(xVal[:, 2:]) == [z1, z2], axis=1))[0] xValInd = xVal[ind][:, :2] yValInd = yVal[ind] valInd = ot.MetaModelValidation(yValInd, resInd.getMetaModel()(xValInd)) error.add(valInd.getResidualSample()) rmseInd = error.computeStandardDeviation()[0] rmseIndList.append(rmseInd) plt.figure() plt.boxplot([rmseLVList, rmseIndList]) plt.xticks([1, 2], ["Mixed continuous-categorical GP", "Independent GPs"]) plt.ylabel("RMSE") .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_categorical_002.png :alt: plot kriging categorical :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_categorical_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(33.722222222222214, 0.5, 'RMSE') .. GENERATED FROM PYTHON SOURCE LINES 369-372 The obtained results show, for this test-case, a better modeling performance when modeling the function as a mixed categorical/continuous function, rather than relying on multiple purely continuous Gaussian processes. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 11.885 seconds) .. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_plot_kriging_categorical.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_kriging_categorical.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_kriging_categorical.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_kriging_categorical.zip `