Note
Go to the end to download the full example code.
Benchmark the Ishigami test function¶
import openturns as ot
import otbenchmark as otb
import openturns.viewer as otv
problem = otb.IshigamiSensitivity()
print(problem)
name = Ishigami
distribution = ComposedDistribution(Uniform(a = -3.14159, b = 3.14159), Uniform(a = -3.14159, b = 3.14159), Uniform(a = -3.14159, b = 3.14159), IndependentCopula(dimension = 3))
function = ParametricEvaluation([X1,X2,X3,a,b]->[sin(X1) + a * sin(X2)^2 + b * X3^4 * sin(X1)], parameters positions=[3,4], parameters=[a : 7, b : 0.1], input positions=[0,1,2])
firstOrderIndices = [0.313905,0.442411,0]
totalOrderIndices = [0.557589,0.442411,0.243684]
distribution = problem.getInputDistribution()
model = problem.getFunction()
Exact first and total order
exact_first_order = problem.getFirstOrderIndices()
exact_first_order
exact_total_order = problem.getTotalOrderIndices()
exact_total_order
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()
full_sample = ot.Sample(size, 1 + dimension)
full_sample[:, range(dimension)] = inputDesign
full_sample[:, dimension] = outputDesign
full_description = list(inputDesign.getDescription())
full_description.append(outputDesign.getDescription()[0])
full_sample.setDescription(full_description)
marginal_distribution = ot.ComposedDistribution(
[
ot.KernelSmoothing().build(full_sample.getMarginal(i))
for i in range(1 + dimension)
]
)
clouds = ot.VisualTest.DrawPairsMarginals(full_sample, marginal_distribution)
_ = otv.View(clouds, figure_kw={"figsize": (6.0, 6.0)})
output_distribution = ot.KernelSmoothing().build(outputDesign)
_ = otv.View(output_distribution.drawPDF())
Perform sensitivity analysis¶
Create X/Y data
ot.RandomGenerator.SetSeed(0)
size = 10000
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 : 10000
Computed first order = [0.302745,0.460846,0.0066916]
Exact first order = [0.313905,0.442411,0]
Computed total order = [0.574996,0.427126,0.256689]
Exact total order = [0.557589,0.442411,0.243684]
_ = otv.View(sensitivityAnalysis.draw())
otv.View.ShowAll()
Total running time of the script: (0 minutes 1.037 seconds)