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
|
"""
Use the ANCOVA indices
======================
"""
# %%
# In this example we are going to use the ANCOVA decomposition to estimate
# sensitivity indices from a model with correlated inputs.
#
# ANCOVA allows one to estimate the Sobol' indices, and thanks to a functional
# decomposition of the model it allows one to separate the part of variance
# explained by a variable itself from the part of variance explained by a
# correlation which is due to its correlation with the other input parameters.
#
#
# In theory, ANCOVA indices range is :math:`\left[0; 1\right]` ; the closer to 1 the index is,
# the greater the model response sensitivity to the variable is.
# These indices are a sum of a physical part :math:`S_i^U` and correlated part :math:`S_i^C`.
# The correlation has a weak influence on the contribution of :math:`X_i`, if :math:`|S_i^C|`
# is low and :math:`S_i` is close to :math:`S_i^U`.
# On the contrary, the correlation has a strong influence on the contribution of
# the input :math:`X_i`, if :math:`|S_i^C|` is high and :math:`S_i` is close to :math:`S_i^C`.
#
# The ANCOVA indices :math:`S_i` evaluate the importance of one variable at a time
# (:math:`d` indices stored, with :math:`d` the input dimension of the model).
# The :math:`d` uncorrelated parts of variance of the output due to each input :math:`S_i^U`
# and the effects of the correlation are represented by the indices resulting
# from the subtraction of these two previous lists.
#
# To evaluate the indices, the library needs of a functional chaos result
# approximating the model response with uncorrelated inputs and a sample with
# correlated inputs used to compute the real values of the output.
#
# %%
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)
# %%
# Create the model (x1,x2) --> (y) = (4.*x1+5.*x2)
model = ot.SymbolicFunction(["x1", "x2"], ["4.*x1+5.*x2"])
# %%
# Create the input independent joint distribution
distribution = ot.Normal(2)
distribution.setDescription(["X1", "X2"])
# %%
# Create the correlated input distribution
S = ot.CorrelationMatrix(2)
S[1, 0] = 0.3
R = ot.NormalCopula.GetCorrelationFromSpearmanCorrelation(S)
copula = ot.NormalCopula(R)
distribution_corr = ot.JointDistribution([ot.Normal()] * 2, copula)
# %%
# ANCOVA needs a functional decomposition of the model
enumerateFunction = ot.LinearEnumerateFunction(2)
productBasis = ot.OrthogonalProductPolynomialFactory(
[ot.HermiteFactory()] * 2, enumerateFunction
)
adaptiveStrategy = ot.FixedStrategy(
productBasis, enumerateFunction.getStrataCumulatedCardinal(4)
)
samplingSize = 250
experiment = ot.MonteCarloExperiment(distribution, samplingSize)
X = experiment.generate()
Y = model(X)
algo = ot.FunctionalChaosAlgorithm(X, Y, distribution, adaptiveStrategy)
algo.run()
result = algo.getResult()
# %%
# Create the input sample taking account the correlation
size = 2000
sample = distribution_corr.getSample(size)
# %%
# Perform the decomposition
ancova = ot.ANCOVA(result, sample)
# Compute the ANCOVA indices (first order and uncorrelated indices are computed together)
indices = ancova.getIndices()
# Retrieve uncorrelated indices
uncorrelatedIndices = ancova.getUncorrelatedIndices()
# Retrieve correlated indices:
correlatedIndices = indices - uncorrelatedIndices
# %%
# Print Sobol' indices, the physical part, and the correlation part
print("ANCOVA indices ", indices)
print("ANCOVA uncorrelated indices ", uncorrelatedIndices)
print("ANCOVA correlated indices ", correlatedIndices)
# %%
graph = ot.SobolIndicesAlgorithm.DrawImportanceFactors(
indices, distribution.getDescription(), "ANCOVA indices (Sobol')"
)
view = viewer.View(graph)
# %%
graph = ot.SobolIndicesAlgorithm.DrawImportanceFactors(
uncorrelatedIndices,
distribution.getDescription(),
"ANCOVA uncorrelated indices\n(part of physical variance in the model)",
)
view = viewer.View(graph)
# %%
graph = ot.SobolIndicesAlgorithm.DrawImportanceFactors(
correlatedIndices,
distribution.getDescription(),
"ANCOVA correlated indices\n(part of variance due to the correlation)",
)
view = viewer.View(graph)
plt.show()
|
