Create a polynomial chaos metamodel from a data set

In this example, we create a polynomial chaos expansion (PCE) using a data set. More precisely, given a pair of input and output samples, we create a PCE without the knowledge of the input distribution. In this example, we use a relatively small sample size equal to 80.

In this example we create a global approximation of a model response using polynomial chaos expansion.

Let be the function defined by:

for any .

We assume that

and that and are independent.

An interesting point in this example is that the output is multivariate. This is why we are going to use the getMarginal method in the script in order to select the output marginal that we want to manage.

Simulate input and output samples

import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt

ot.Log.Show(ot.Log.NONE)

We first create the function model.

ot.RandomGenerator.SetSeed(2)
dimension = 2
input_names = ["x1", "x2"]
formulas = ["cos(x1 + x2)", "(x2 + 1) * exp(x1)"]
model = ot.SymbolicFunction(input_names, formulas)

Then we create a sample inputSample and compute the corresponding output sample outputSample.

distribution = ot.Normal(dimension)
samplesize = 80
inputSample = distribution.getSample(samplesize)
outputSample = model(inputSample)

Create the PCE

Create a functional chaos model. The algorithm needs to fit a distribution on the input sample. To do this, the algorithm in FunctionalChaosAlgorithm uses the BuildDistribution static method to fit the distribution to the input sample. Please read Fit a distribution from an input sample for an example of this method. The algorithm does this automatically using the Lilliefors test. In order to make the algorithm a little faster, we reduce the value of the sample size used in the Lilliefors test.

ot.ResourceMap.SetAsUnsignedInteger("FittingTest-LillieforsMaximumSamplingSize", 50)

The main topic of this example is to introduce the next constructor of FunctionalChaosAlgorithm. Notice that the only input arguments are the input and output sample.

algo = ot.FunctionalChaosAlgorithm(inputSample, outputSample)
algo.run()
result = algo.getResult()
result

FunctionalChaosResult

• input dimension: 2
• output dimension: 2
• distribution dimension: 2
• transformation: 2 -> 2
• inverse transformation: 2 -> 2
• orthogonal basis dimension: 2
• indices size: 33
• relative errors: [0.000423845,3.97704e-08]
• residuals: [0.00156354,0.000147243]

Index	Multi-index	Coeff.#0	Coeff.#1
0	[0,0]	-0.9778362	2.494928
1	[1,0]	3.946319	1.310505
2	[0,1]	-0.5668731	2.058245
3	[2,0]	-6.595425	3.31883
4	[0,2]	0.05754212	0.1086024
5	[3,0]	7.262131	-0.841792
6	[0,3]	-0.9758324	-0.12278
7	[1,1]	-0.7887584	2.418963
8	[4,0]	-6.489266	1.935301
9	[0,4]	0.2844189	0.1285879
10	[5,0]	5.352722	-1.168838
11	[0,5]	-0.8491851	-0.1064447
12	[2,1]	0.3063163	1.643328
13	[1,2]	-0.0281957	0.01273302
14	[6,0]	-3.800862	0.9235032
15	[0,6]	0.004040826	0.1013839
16	[7,0]	2.402248	-0.5406551
17	[0,7]	-0.441342	-0.0634571
18	[3,1]	-0.003325215	1.022703
19	[1,3]	-0.2267971	0.002919773
20	[2,2]	0.4659559	-0.001704627
21	[8,0]	-1.272689	0.2902903
22	[0,8]	-0.03341318	0.05202352
23	[4,1]	0.1240007	0.2562661
24	[1,4]	-0.04681154	-0.006325388
25	[9,0]	0.5059943	-0.1094002
26	[0,9]	-0.1155697	-0.01848551
27	[3,2]	-0.04218489	0.02369182
28	[2,3]	-0.05014101	-0.01482048
29	[10,0]	-0.1825852	0.03945734
30	[0,10]	-0.01467434	0.01355966
31	[5,1]	-0.08889085	0.1317192
32	[1,5]	-0.2461735	0.001139069

Not all coefficients are selected in this PCE. Indeed, the default constructor of FunctionalChaosAlgorithm creates a sparse PCE. Please read Create a full or sparse polynomial chaos expansion for more details on this topic.

Get the metamodel.

metamodel = result.getMetaModel()

Plot the second output of our model depending on with . In order to do this, we create a ParametricFunction and set the value of . Then we use the getMarginal method to extract the second output (which index is equal to 1).

x1index = 0
x1value = 0.5
x2min = -3.0
x2max = 3.0
outputIndex = 1
metamodelParametric = ot.ParametricFunction(metamodel, [x1index], [x1value])
graph = metamodelParametric.getMarginal(outputIndex).draw(x2min, x2max)
graph.setLegends(["Metamodel"])
modelParametric = ot.ParametricFunction(model, [x1index], [x1value])
curve = modelParametric.getMarginal(outputIndex).draw(x2min, x2max).getDrawable(0)
curve.setColor("red")
curve.setLegend("Model")
graph.add(curve)
graph.setLegendPosition("lower right")
graph.setXTitle("X2")
graph.setTitle("Metamodel Validation, output #%d" % (outputIndex))
view = viewer.View(graph)

We see that the metamodel fits approximately to the model, except perhaps for extreme values of . However, there is a better way of globally validating the metamodel, using the MetaModelValidation on a validation design of experiment.

n_valid = 100
inputTest = distribution.getSample(n_valid)
outputTest = model(inputTest)

Plot the corresponding validation graphics.

val = ot.MetaModelValidation(outputTest, metamodel(inputTest))
R2 = val.computeR2Score()
graph = val.drawValidation()
graph.setTitle("Metamodel validation R2=" + str(R2))
view = viewer.View(graph)

The coefficient of predictivity is not extremely satisfactory for the first output, but is would be sufficient for a central dispersion study. The second output has a much more satisfactory R2: only one single extreme point is far from the diagonal of the graphics.

Compute and print Sobol’ indices

chaosSI = ot.FunctionalChaosSobolIndices(result)
chaosSI

FunctionalChaosSobolIndices

• input dimension: 2
• output dimension: 2
• basis size: 33
• mean: [-0.977836,2.49493]
• std-dev: [14.4244,5.81255]

Output marginal: 0

Input	Variable	Sobol' index	Total index
0	X0	0.983831	0.989001
1	X1	0.010999	0.016169

Index	Multi-index	Part of variance
5	[3,0]	0.253472
3	[2,0]	0.209068
8	[4,0]	0.202392
10	[5,0]	0.137706
1	[1,0]	0.074849
14	[6,0]	0.069433
16	[7,0]	0.027736

Output marginal: 1

Input	Variable	Sobol' index	Total index
0	X0	0.585905	0.872471
1	X1	0.127529	0.414095

Index	Multi-index	Part of variance
3	[2,0]	0.326015
7	[1,1]	0.173191
2	[0,1]	0.125390
8	[4,0]	0.110857
12	[2,1]	0.079931
1	[1,0]	0.050833
10	[5,0]	0.040437
18	[3,1]	0.030958
14	[6,0]	0.025243
5	[3,0]	0.020974

Let us analyse the results of this global sensitivity analysis.

• We see that the first output involves significant multi-indices with higher marginal degree.

• For the second output, the first variable is especially significant, with a significant contribution of the interactions. The contribution of the interactions are very significant in this model.

Draw Sobol’ indices.

sensitivityAnalysis = ot.FunctionalChaosSobolIndices(result)
first_order = [sensitivityAnalysis.getSobolIndex(i) for i in range(dimension)]
total_order = [sensitivityAnalysis.getSobolTotalIndex(i) for i in range(dimension)]

input_names = model.getInputDescription()
graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(input_names, first_order, total_order)
graph.setLegendPosition("center")
view = viewer.View(graph)

Testing the sensitivity to the degree

With the specific constructor of FunctionalChaosAlgorithm that we use, the FunctionalChaosAlgorithm-MaximumTotalDegree in the ResourceMap configure the maximum degree explored by the algorithm. This degree is a trade-off.

• If the maximum degree is too low, the polynomial may miss some coefficients so that the quality is lower than possible.

• If the maximum degree is too large, the number of coefficients to explore is too large, so that the coefficients might be poorly estimated.

This is why the following for loop explores various degrees to see the sensitivity of the metamodel predictivity depending on the degree.

The default value of this parameter is 10.

ot.ResourceMap.GetAsUnsignedInteger("FunctionalChaosAlgorithm-MaximumTotalDegree")

10

This is why we explore the values from 1 to 15.

degrees = range(1, 12)
r2 = ot.Sample(len(degrees), 2)
for maximumDegree in degrees:
    ot.ResourceMap.SetAsUnsignedInteger(
        "FunctionalChaosAlgorithm-MaximumTotalDegree", maximumDegree
    )
    print("Maximum total degree =", maximumDegree)
    algo = ot.FunctionalChaosAlgorithm(inputSample, outputSample)
    algo.run()
    result = algo.getResult()
    metamodel = result.getMetaModel()
    for outputIndex in range(2):
        val = ot.MetaModelValidation(
            outputTest[:, outputIndex], metamodel.getMarginal(outputIndex)(inputTest)
        )
        r2Value = min(1.0, max(0.0, val.computeR2Score()[0]))  # Get lucky.
        r2[maximumDegree - degrees[0], outputIndex] = r2Value

Maximum total degree = 1
Maximum total degree = 2
Maximum total degree = 3
Maximum total degree = 4
Maximum total degree = 5
Maximum total degree = 6
Maximum total degree = 7
Maximum total degree = 8
Maximum total degree = 9
Maximum total degree = 10
Maximum total degree = 11

graph = ot.Graph("Predictivity", "Total degree", "R2", True)
cloud = ot.Cloud([[d] for d in degrees], r2[:, 0])
cloud.setLegend("Output #0")
cloud.setPointStyle("bullet")
graph.add(cloud)
cloud = ot.Cloud([[d] for d in degrees], r2[:, 1])
cloud.setLegend("Output #1")
cloud.setColor("red")
cloud.setPointStyle("bullet")
graph.add(cloud)
graph.setLegendPosition("upper right")
view = viewer.View(graph, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"})
plt.subplots_adjust(right=0.7)

plt.show()

We see that the R2 score increases then gets constant or decreases. A low total polynomial degree is not sufficient to describe the first output with good predictivity. However, the coefficient of predictivity can decrease when the total degree gets greater. The predictivity of the second output seems to be much less satisfactory: a little more work would be required to improve the metamodel.

In this situation, the following methods may be used.

• Since the distribution of the input is known, we may want to give this information to the FunctionalChaosAlgorithm. This prevents the algorithm from trying to fit the input distribution which best fit to the data.

• We may want to customize the adaptiveStrategy by selecting an enumerate function which best fit to this particular example. In this specific example, however, the interactions plays a great role so that the linear enumerate function may provide better results than the hyperbolic rule.

• We may want to customize the projectionStrategy by selecting a method to compute the coefficient which improves the estimation. For example, it might be interesting to try an integration rule instead of the least squares method. Notice that a specific design of experiment is required in this case.

Reset default settings

ot.ResourceMap.Reload()