Benchmark the Morris test function

import openturns as ot
import otbenchmark as otb
import openturns.viewer as otv
problem = otb.MorrisSensitivity()
print(problem)
name = Morris
distribution = ComposedDistribution(Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), Uniform(a = 0, b = 1), IndependentCopula(dimension = 20))
function = class=PythonEvaluation name=MorrisFunction
firstOrderIndices = [0.08,0.08,0.06,0.08,0.06,0.13,0.06,0.13,0.13,0.12,0,0,0,0,0,0,0,0,0,0]#20
totalOrderIndices = [0.11,0.11,0.06,0.11,0.06,0.13,0.06,0.13,0.13,0.12,0,0,0,0,0,0,0,0,0,0]#20
distribution = problem.getInputDistribution()
model = problem.getFunction()

Exact first and total order

exact_first_order = problem.getFirstOrderIndices()
exact_first_order
class=Point name=Unnamed dimension=20 values=[0.08,0.08,0.06,0.08,0.06,0.13,0.06,0.13,0.13,0.12,0,0,0,0,0,0,0,0,0,0]


exact_total_order = problem.getTotalOrderIndices()
exact_total_order
class=Point name=Unnamed dimension=20 values=[0.11,0.11,0.06,0.11,0.06,0.13,0.06,0.13,0.13,0.12,0,0,0,0,0,0,0,0,0,0]


Plot the function

Create X/Y data

ot.RandomGenerator.SetSeed(0)
size = 200
inputDesign = ot.MonteCarloExperiment(distribution, size).generate()
outputDesign = model(inputDesign)
dimension = distribution.getDimension()
nbcolumns = 4
nbrows = int(dimension / nbcolumns)
grid = ot.GridLayout(nbrows, nbcolumns)
inputDescription = distribution.getDescription()
outputDescription = model.getOutputDescription()[0]
index = 0
for i in range(nbrows):
    for j in range(nbcolumns):
        graph = ot.Graph(
            "n=%d" % (size), inputDescription[index], outputDescription, True, ""
        )
        sample = ot.Sample(size, 2)
        sample[:, 0] = inputDesign[:, index]
        sample[:, 1] = outputDesign[:, 0]
        cloud = ot.Cloud(sample)
        graph.add(cloud)
        grid.setGraph(i, j, graph)
        index += 1
_ = otv.View(grid, figure_kw={"figsize": (10.0, 10.0)})

# %
output_distribution = ot.KernelSmoothing().build(outputDesign)
_ = otv.View(output_distribution.drawPDF())
  • , n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200, n=200
  • plot morris sensitivity

Perform sensitivity analysis

Create X/Y data

ot.RandomGenerator.SetSeed(0)
size = 30
inputDesign = ot.SobolIndicesExperiment(distribution, size).generate()
outputDesign = model(inputDesign)

Compute first order indices using the Saltelli estimator

sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(inputDesign, outputDesign, size)
computed_first_order = sensitivityAnalysis.getFirstOrderIndices()
computed_total_order = sensitivityAnalysis.getTotalOrderIndices()

Compare with exact results

print("Sample size : ", size)
# First order
# Compute absolute error (the LRE cannot be computed,
# because S can be zero)
print("Computed first order = ", computed_first_order)
print("Exact first order = ", exact_first_order)
# Total order
print("Computed total order = ", computed_total_order)
print("Exact total order = ", exact_total_order)
Sample size :  30
Computed first order =  [-0.133001,-0.284659,-0.244261,-0.013848,0.0717474,-0.189384,-0.192788,-0.232371,0.176264,-0.119774,-0.220786,-0.275918,-0.232678,-0.25769,-0.276469,-0.278563,-0.261348,-0.287263,-0.242783,-0.248379]#20
Exact first order =  [0.08,0.08,0.06,0.08,0.06,0.13,0.06,0.13,0.13,0.12,0,0,0,0,0,0,0,0,0,0]#20
Computed total order =  [0.14052,0.0954675,-0.0522002,0.153156,0.0131494,0.113636,0.109581,0.150205,0.169807,0.268959,0.0473076,0.00448006,0.00150032,0.00855239,-0.0238276,-0.0252896,0.0187677,-0.00886828,-0.0161144,-0.00844209]#20
Exact total order =  [0.11,0.11,0.06,0.11,0.06,0.13,0.06,0.13,0.13,0.12,0,0,0,0,0,0,0,0,0,0]#20
_ = otv.View(sensitivityAnalysis.draw())
Sobol' indices - SaltelliSensitivityAlgorithm
otv.View.ShowAll()

Total running time of the script: (0 minutes 14.886 seconds)