1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
|
#! /usr/bin/env python
import openturns as ot
from openturns.usecases import ishigami_function
ot.TESTPREAMBLE()
input_dimension = 3
output_dimension = 1
formula = [
"sin(pi_*X1)+7*sin(pi_*X2)*sin(pi_*X2)+0.1*((pi_*X3)*(pi_*X3)*(pi_*X3)*(pi_*X3))*sin(pi_*X1)"
]
model = ot.SymbolicFunction(["X1", "X2", "X3"], formula)
distribution = ot.JointDistribution([ot.Uniform(-1.0, 1.0)] * input_dimension)
# Size of simulation
size = 10000
# Test with the various implementation methods
methods = ["Saltelli", "Jansen", "MauntzKucherenko", "Martinez"]
# Test the different sampling options
samplings = ["MonteCarlo", "LHS", "QMC"]
# Generate input/output designs
computeSO = True
# Case 1 : Estimation of sensitivity using estimator and no bootstrap
for method in methods:
for sampling in samplings:
ot.ResourceMap.SetAsString("SobolIndicesExperiment-SamplingMethod", sampling)
sensitivity_algorithm = eval(
"ot."
+ method
+ "SensitivityAlgorithm(distribution, size, model, computeSO)"
)
print("Method of evaluation=", method)
print("Method of sampling=", sampling)
# Get first order indices
fo = sensitivity_algorithm.getFirstOrderIndices()
print("First order indices = ", fo)
# Get total order indices
to = sensitivity_algorithm.getTotalOrderIndices()
print("Total order indices = ", to)
# Get the confidence interval thanks to Bootstrap
nr_bootstrap = 100
confidence_level = 0.95
sensitivity_algorithm.setBootstrapSize(nr_bootstrap)
sensitivity_algorithm.setConfidenceLevel(confidence_level)
sensitivity_algorithm.setUseAsymptoticDistribution(False)
interval_fo = sensitivity_algorithm.getFirstOrderIndicesInterval()
interval_to = sensitivity_algorithm.getTotalOrderIndicesInterval()
print("bootstrap intervals")
print("First order indices interval = ", interval_fo)
print("Total order indices interval = ", interval_to)
# Asymptotic confidence interval
sensitivity_algorithm.setUseAsymptoticDistribution(True)
interval_fo_asymptotic = sensitivity_algorithm.getFirstOrderIndicesInterval()
interval_to_asymptotic = sensitivity_algorithm.getTotalOrderIndicesInterval()
print("asymptotic intervals:")
print(
"First order indices distribution = ",
sensitivity_algorithm.getFirstOrderIndicesDistribution(),
)
print(
"Total order indices distribution = ",
sensitivity_algorithm.getTotalOrderIndicesDistribution(),
)
print("First order indices interval = ", interval_fo_asymptotic)
print("Total order indices interval = ", interval_to_asymptotic)
# with experiment
sequence = ot.SobolSequence(input_dimension)
experiment = ot.LowDiscrepancyExperiment(
sequence, ot.JointDistribution([ot.Uniform(0.0, 1.0)] * input_dimension), size
)
sensitivity_algorithm = ot.SaltelliSensitivityAlgorithm(experiment, model)
print(sensitivity_algorithm.getFirstOrderIndices())
# multi variate model
model_aggregated = ot.SymbolicFunction(
["X1", "X2", "X3"],
["2*X1 + X2 - 3*X3 + 0.3*X1*X2", "-5*X1 + 4*X2 - 0.8*X2*X3 + 2*X3"],
)
distribution_aggregated = ot.JointDistribution([ot.Uniform()] * 3)
inputDesign = ot.SobolIndicesExperiment(distribution_aggregated, size).generate()
outputDesign = model_aggregated(inputDesign)
# Case 1 : Estimation of sensitivity using estimator and no bootstrap
for method in methods:
sensitivity_algorithm = eval(
"ot." + method + "SensitivityAlgorithm(inputDesign, outputDesign, size)"
)
print("Method of evaluation=", method)
# Get first order indices
fo = sensitivity_algorithm.getAggregatedFirstOrderIndices()
print("Aggregated first order indices = ", fo)
# Get total order indices
to = sensitivity_algorithm.getAggregatedTotalOrderIndices()
print("Aggregated total order indices = ", to)
# Get the confidence interval thanks to Bootstrap
nr_bootstrap = 100
confidence_level = 0.95
# sensitivity_algorithm = ot.MartinezSensitivityAlgorithm(
# inputDesign, outputDesign, size)
sensitivity_algorithm.setBootstrapSize(nr_bootstrap)
sensitivity_algorithm.setConfidenceLevel(confidence_level)
sensitivity_algorithm.setUseAsymptoticDistribution(False)
interval_fo = sensitivity_algorithm.getFirstOrderIndicesInterval()
interval_to = sensitivity_algorithm.getTotalOrderIndicesInterval()
print("bootstrap intervals")
print("Aggregated first order indices interval = ", interval_fo)
print("Aggregated total order indices interval = ", interval_to)
# Asymptotic confidence interval
sensitivity_algorithm.setUseAsymptoticDistribution(True)
interval_fo_asymptotic = sensitivity_algorithm.getFirstOrderIndicesInterval()
interval_to_asymptotic = sensitivity_algorithm.getTotalOrderIndicesInterval()
print("asymptotic intervals:")
print("Aggregated first order indices interval = ", interval_fo_asymptotic)
print("Aggregated total order indices interval = ", interval_to_asymptotic)
# special case in dim=2
ot.ResourceMap.SetAsString("SobolIndicesExperiment-SamplingMethod", "MonteCarlo")
ot.RandomGenerator.SetSeed(0)
distribution = ot.JointDistribution([ot.Uniform()] * 2)
size = 1000
model = ot.SymbolicFunction(["X1", "X2"], ["2*X1 + X2"])
sensitivity_algorithm = ot.SaltelliSensitivityAlgorithm(distribution, size, model, True)
print(sensitivity_algorithm.getSecondOrderIndices())
ot.RandomGenerator.SetSeed(0)
experiment = ot.MonteCarloExperiment(distribution, size)
sensitivity_algorithm = ot.SaltelliSensitivityAlgorithm(experiment, model, True)
print(sensitivity_algorithm.getSecondOrderIndices())
ot.RandomGenerator.SetSeed(0)
x = ot.SobolIndicesExperiment(distribution, size, True).generate()
y = model(x)
sensitivity_algorithm = ot.SaltelliSensitivityAlgorithm(x, y, size)
print(sensitivity_algorithm.getSecondOrderIndices())
# null contribution case: X3 not in output formula
model = ot.SymbolicFunction(["X1", "X2", "X3"], ["10+3*X1+X2"])
distribution = ot.JointDistribution([ot.Uniform(-1.0, 1.0)] * input_dimension)
size = 10000
for method in methods:
sensitivity_algorithm = eval(
"ot." + method + "SensitivityAlgorithm(distribution, size, model, False)"
)
sensitivity_algorithm.setUseAsymptoticDistribution(True)
fo = sensitivity_algorithm.getFirstOrderIndices()
print("First order indices = ", fo)
to = sensitivity_algorithm.getTotalOrderIndices()
# print("Total order indices = ", to)
interval_fo = sensitivity_algorithm.getFirstOrderIndicesInterval()
interval_to = sensitivity_algorithm.getTotalOrderIndicesInterval()
print("Aggregated first order indices interval = ", interval_fo)
# print("Aggregated total order indices interval = ", interval_to)
# setDesign must reset results across runs
ot.Log.Show(ot.Log.NONE)
im = ishigami_function.IshigamiModel()
exact_first_order = ot.Point([im.S1, im.S2, im.S3])
exact_total_order = ot.Point([im.ST1, im.ST2, im.ST3])
sobolIndicesAlgorithmB = ot.SaltelliSensitivityAlgorithm()
for sample_size in [100, 1000, 10000]:
print("Size:", sample_size)
# Method A : classical
X = ot.SobolIndicesExperiment(im.distribution, sample_size).generate()
Y = im.model(X)
sobolIndicesAlgorithmA = ot.SaltelliSensitivityAlgorithm(X, Y, sample_size)
computed_first_orderA = sobolIndicesAlgorithmA.getFirstOrderIndices()
computed_total_orderA = sobolIndicesAlgorithmA.getTotalOrderIndices()
first_error = computed_first_orderA - exact_first_order
total_error = computed_total_orderA - exact_total_order
print("Method A :", first_error, total_error)
# Method B : setDesign
sobolIndicesAlgorithmB.setDesign(X, Y, sample_size)
computed_first_orderB = sobolIndicesAlgorithmB.getFirstOrderIndices()
computed_total_orderB = sobolIndicesAlgorithmB.getTotalOrderIndices()
first_error = computed_first_orderB - exact_first_order
total_error = computed_total_orderB - exact_total_order
print("Method B :", first_error, total_error)
assert computed_first_orderB == computed_first_orderA, "wrong first"
assert computed_total_orderB == computed_total_orderA, "wrong total"
|
