Polynomial chaos exploitation¶

In this example we are going to create a global approximation of a model response using functional chaos and expose the associated results:

the composed model: $h$ : $\underline{Z}^{\strut} \longrightarrow \underline{Y} = g \circ T^{-1}(\underline{Z})$ , which is the model of the reduced variables $\underline{Z}$ . We have $\displaystyle h = \sum_{k \in \mathbb N} \underline{\alpha}_k \Psi_k$ ,
the coefficients of the polynomial approximation : $(\underline{\alpha}_k)_{k \in K}$ ,
the composed meta model: $\hat{h}$ , which is the model of the reduced variables reduced to the truncated multivariate basis $(\Psi_k)_{k \in K}$ . We have $\displaystyle \hat{h} = \sum_{k \in K} \underline{\alpha}_k \Psi_k$ ,
the meta model: $\displaystyle \hat{g} : \underline{X} \longrightarrow Y = \hat{h} \circ T(\underline{X})$ which is the polynomial chaos approximation as a Function. We have $\displaystyle \hat{g} = \sum_{k \in K} \underline{\alpha}_k \Psi_k \circ T$ ,
the truncated multivariate basis : $(\Psi_k)_{k \in K}$ ,
the indices $K$ ,
the composition of each polynomial of the truncated multivariate basis $\Psi_k$ ,
the distribution $\mu$ of the transformed variables $\underline{Z}$ ,

[1]:

from __future__ import print_function
import openturns as ot

[2]:

# prepare some X/Y data
ot.RandomGenerator.SetSeed(0)
dimension = 2
input_names = ['x1', 'x2']
formulas = ['cos(x1 + x2)', '(x2 + 1) * exp(x1 - 2 * x2)']
model = ot.SymbolicFunction(input_names, formulas)
distribution = ot.Normal(dimension)
x  = distribution.getSample(30)
y = model(x)

[3]:

# create a functional chaos algorithm
algo = ot.FunctionalChaosAlgorithm(x, y)
algo.run()

[4]:

# Stream out the result
result = algo.getResult()

[5]:

# Get the polynomial chaos coefficients
# All the coefficients
result.getCoefficients()

[5]:

	v0	v1
0	0.42562857011464145	0.19814108237713063
1	0.0	0.470485758115012
2	-0.09932089428332584	0.7885721936793412
3	-0.24625085236020092	0.3232951290142421
4	-0.19204843536289498	-1.8462421285361612
5	0.0	1.4246957691984192
6	-0.47643611347170484	0.0
7	0.0	-1.029321781268946
8	0.0	-0.1355602355790481
9	0.0	-0.28170119639243424
10	-0.0249183372934895	0.0
11	0.0	-0.30064668963377356
12	0.07558242403322926	0.0
13	0.0	0.34014180081319967
14	0.18874495042818668	0.0
15	0.0	0.03466594712788605

[6]:

# The coefficients of marginal i
i = 1
result.getCoefficients()[i]

[6]:

[0,0.470486]

[7]:

# Get the indices of the selected polynomials : K
subsetK = result.getIndices()
subsetK

[7]:

[0,1,2,3,4,6,7,9,12,13,16,17,18,19,20,21]#16

[8]:

# Get the composition of the polynomials
# of the truncated multivariate basis
for i in range(subsetK.getSize()):
    print("Polynomial number ", i, " in truncated basis <-> polynomial number ",
          subsetK[i], " = ", ot.LinearEnumerateFunction(dimension)(subsetK[i]), " in complete basis")

Polynomial number  0  in truncated basis <-> polynomial number  0  =  [0,0]  in complete basis
Polynomial number  1  in truncated basis <-> polynomial number  1  =  [1,0]  in complete basis
Polynomial number  2  in truncated basis <-> polynomial number  2  =  [0,1]  in complete basis
Polynomial number  3  in truncated basis <-> polynomial number  3  =  [2,0]  in complete basis
Polynomial number  4  in truncated basis <-> polynomial number  4  =  [1,1]  in complete basis
Polynomial number  5  in truncated basis <-> polynomial number  6  =  [3,0]  in complete basis
Polynomial number  6  in truncated basis <-> polynomial number  7  =  [2,1]  in complete basis
Polynomial number  7  in truncated basis <-> polynomial number  9  =  [0,3]  in complete basis
Polynomial number  8  in truncated basis <-> polynomial number  12  =  [2,2]  in complete basis
Polynomial number  9  in truncated basis <-> polynomial number  13  =  [1,3]  in complete basis
Polynomial number  10  in truncated basis <-> polynomial number  16  =  [4,1]  in complete basis
Polynomial number  11  in truncated basis <-> polynomial number  17  =  [3,2]  in complete basis
Polynomial number  12  in truncated basis <-> polynomial number  18  =  [2,3]  in complete basis
Polynomial number  13  in truncated basis <-> polynomial number  19  =  [1,4]  in complete basis
Polynomial number  14  in truncated basis <-> polynomial number  20  =  [0,5]  in complete basis
Polynomial number  15  in truncated basis <-> polynomial number  21  =  [6,0]  in complete basis

[9]:

# Get the multivariate basis
# as a collection of Function
reduced = result.getReducedBasis()

[10]:

# Get the orthogonal basis
orthgBasis = result.getOrthogonalBasis()
orthgBasis

[10]:

class=OrthogonalProductPolynomialFactory univariate polynomial collection=[class=OrthogonalUniVariatePolynomialFamily implementation=class=StandardDistributionPolynomialFactory hasSpecificFamily=true specificFamily=class=OrthogonalUniVariatePolynomialFamily implementation=class=HistogramPolynomialFactory measure=class=Histogram name=Histogram dimension=1 first=-1 width=class=Point name=Unnamed dimension=5 values=[0.4,0.4,0.4,0.4,0.4] height=class=Point name=Unnamed dimension=5 values=[0.25,0.25,0.666667,1.08333,0.25] legendre=class=LegendreFactory measure=class=Uniform name=Uniform dimension=1 a=-1 b=1,class=OrthogonalUniVariatePolynomialFamily implementation=class=StandardDistributionPolynomialFactory hasSpecificFamily=false orthonormalization algorithm=class=OrthonormalizationAlgorithm implementation=class=AdaptiveStieltjesAlgorithm measure=class=Trapezoidal name=Trapezoidal dimension=1 a=-1 b=-0.530738 c=0.0480939 d=1 h=0.775545 monicRecurrenceCoefficients=[class=Point name=Unnamed dimension=3 values=[1,0.098496,0],class=Point name=Unnamed dimension=3 values=[1,0.0055215,-0.184989],class=Point name=Unnamed dimension=3 values=[1,-0.0074096,-0.222231],class=Point name=Unnamed dimension=3 values=[1,0.0105003,-0.237692],class=Point name=Unnamed dimension=3 values=[1,-0.00366255,-0.243879],class=Point name=Unnamed dimension=3 values=[1,0.000574247,-0.242444],class=Point name=Unnamed dimension=3 values=[1,0.000785537,-0.246691],class=Point name=Unnamed dimension=3 values=[1,0.000879275,-0.247369],class=Point name=Unnamed dimension=3 values=[1,-0.00173165,-0.247045]] monicSquaredNorms=class=Point name=Unnamed dimension=9 values=[1,0.184989,0.0411103,0.00977157,0.00238308,0.000577763,0.000142529,3.52572e-05,8.71013e-06] isElliptical=false] measure=class=ComposedDistribution name=ComposedDistribution dimension=2 copula=class=IndependentCopula name=IndependentCopula dimension=2 marginal[0]=class=Histogram name=Histogram dimension=1 first=-1 width=class=Point name=Unnamed dimension=5 values=[0.4,0.4,0.4,0.4,0.4] height=class=Point name=Unnamed dimension=5 values=[0.25,0.25,0.666667,1.08333,0.25] marginal[1]=class=Trapezoidal name=Trapezoidal dimension=1 a=-1 b=-0.530738 c=0.0480939 d=1 h=0.775545

[11]:

# Get the distribution of variables Z
orthgBasis.getMeasure()

[11]:

ComposedDistribution(Histogram(origin = -1, {w0 = 0.4, h0 = 0.25}, {w1 = 0.4, h1 = 0.25}, {w2 = 0.4, h2 = 0.666667}, {w3 = 0.4, h3 = 1.08333}, {w4 = 0.4, h4 = 0.25}), Trapezoidal(a = -1, b = -0.530738, c = 0.0480939, d = 1), IndependentCopula(dimension = 2))

[12]:

# Get the composed model which is the model of the reduced variables Z
result.getComposedModel()

[12]:

(DatabaseEvaluation
input sample :
[ X0 X1 ]
0 : [ 0.608202 -1.26617 ]
1 : [ -0.438266 1.20548 ]
2 : [ -2.18139 0.350042 ]
3 : [ -0.355007 1.43725 ]
4 : [ 0.810668 0.793156 ]
5 : [ -0.470526 0.261018 ]
6 : [ -2.29006 -1.28289 ]
7 : [ -1.31178 -0.0907838 ]
8 : [ 0.995793 -0.139453 ]
9 : [ -0.560206 0.44549 ]
10 : [ 0.322925 0.445785 ]
11 : [ -1.03808 -0.856712 ]
12 : [ 0.473617 -0.125498 ]
13 : [ 0.351418 1.78236 ]
14 : [ 0.0702074 -0.781366 ]
15 : [ -0.721533 -0.241223 ]
16 : [ -1.78796 0.40136 ]
17 : [ 1.36783 1.00434 ]
18 : [ 0.741548 -0.0436123 ]
19 : [ 0.539345 0.29995 ]
20 : [ 0.407717 -0.485112 ]
21 : [ -0.382992 -0.752817 ]
22 : [ 0.257926 1.96876 ]
23 : [ -0.671291 1.85579 ]
24 : [ 0.0521593 0.790446 ]
25 : [ 0.716353 -0.743622 ]
26 : [ 0.184356 -1.53073 ]
27 : [ 0.655027 0.538071 ]
28 : [ 1.73821 -0.958722 ]
29 : [ 0.377922 -0.181004 ]
output sample :
[ y0 y1 ]
0 : [ 0.791234 -6.153 ]
1 : [ 0.719848 0.127674 ]
2 : [ -0.257609 0.075673 ]
3 : [ 0.46935 0.0964592 ]
4 : [ -0.0330217 0.825582 ]
5 : [ 0.978133 0.467366 ]
6 : [ -0.9084 -0.372691 ]
7 : [ 0.167439 0.293644 ]
8 : [ 0.655206 3.07871 ]
9 : [ 0.993427 0.338667 ]
10 : [ 0.718808 0.818737 ]
11 : [ -0.318354 0.28152 ]
12 : [ 0.940016 1.80491 ]
13 : [ -0.533709 0.111917 ]
14 : [ 0.757606 1.11916 ]
15 : [ 0.571259 0.59742 ]
16 : [ 0.183152 0.105058 ]
17 : [ -0.718312 1.05597 ]
18 : [ 0.76617 2.19061 ]
19 : [ 0.667988 1.22357 ]
20 : [ 0.997007 2.04242 ]
21 : [ 0.421399 0.759585 ]
22 : [ -0.609865 0.0749114 ]
23 : [ 0.376759 0.0356671 ]
24 : [ 0.665521 0.388187 ]
25 : [ 0.999628 2.32215 ]
26 : [ 0.222539 -13.6308 ]
27 : [ 0.368781 1.00946 ]
28 : [ 0.711272 1.59716 ]
29 : [ 0.980674 1.71644 ])o(((| y0 = Normal(mu = 0, sigma = 1) -> y0 : Histogram(origin = -2.29006, {w0 = 0.805655, h0 = 0.124123}, {w1 = 0.805655, h1 = 0.124123}, {w2 = 0.805655, h2 = 0.330994}, {w3 = 0.805655, h3 = 0.537865}, {w4 = 0.805655, h4 = 0.124123})
| y1 = Normal(mu = 0, sigma = 1) -> y1 : Trapezoidal(a = -1.75839, b = -0.781359, c = 0.423804, d = 2.40573)
)o(class=LinearEvaluation name=Unnamed center=[0,0] constant=[0,0] linear=[[ 1 0.0108329 ]
[ 0 0.999941 ]]))o(| y0 = Histogram(origin = -1, {w0 = 0.4, h0 = 0.25}, {w1 = 0.4, h1 = 0.25}, {w2 = 0.4, h2 = 0.666667}, {w3 = 0.4, h3 = 1.08333}, {w4 = 0.4, h4 = 0.25}) -> y0 : Normal(mu = 0, sigma = 1)
| y1 = Trapezoidal(a = -1, b = -0.530738, c = 0.0480939, d = 1) -> y1 : Normal(mu = 0, sigma = 1)
))

[13]:

# Get the composed meta model which is the model of the reduced variables Z
# within the reduced polynomials basis
result.getComposedMetaModel()

[13]:

[0.425629,0.198141] + [0,0.470486] * (-0.29173 + 2.18797 * x0) + [-0.0993209,0.788572] * (0.229006 + 2.32502 * x1) + [-0.246251,0.323295] * (-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) + [-0.192048,-1.84624] * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2) + [0,1.4247] * (-0.207607 - 4.12182 * x1 + 0.977306 * x1^2 + 10.1162 * x1^3) + [-0.476436,0] * ((-0.29173 + 2.18797 * x0) * (0.229006 + 2.32502 * x1)) + [0,-1.02932] * (0.893658 - 1.0145 * x1 - 13.1947 * x1^2 + 2.19409 * x1^3 + 20.4848 * x1^4) + [0,-0.13556] * ((-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) * (0.229006 + 2.32502 * x1)) + [0,-0.281701] * ((-0.29173 + 2.18797 * x0) * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2)) + [-0.0249183,0] * (0.284385 - 7.63723 * x0 - 4.67858 * x0^2 + 78.3356 * x0^3 + 9.59256 * x0^4 - 185.818 * x0^5 - 5.17002 * x0^6 + 120.521 * x0^7) + [0,-0.300647] * (-0.202696 - 7.72688 * x1 + 5.82995 * x1^2 + 87.3021 * x1^3 - 20.8178 * x1^4 - 231.97 * x1^5 + 17.6506 * x1^6 + 168.413 * x1^7) + [0.0755824,0] * ((0.400359 - 4.06203 * x0 - 0.288198 * x0^2 + 7.74638 * x0^3) * (0.229006 + 2.32502 * x1)) + [0,0.340142] * ((-0.29173 + 2.18797 * x0) * (-0.207607 - 4.12182 * x1 + 0.977306 * x1^2 + 10.1162 * x1^3)) + [0.188745,0] * ((-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2)) + [0,0.0346659] * (0.95105 + 2.92377 * x0 - 42.0158 * x0^2 - 17.3514 * x0^3 + 250.399 * x0^4 + 29.114 * x0^5 - 454.203 * x0^6 - 14.8305 * x0^7 + 250.574 * x0^8)

[14]:

# Get the meta model which is the composed meta model combined with the
# iso probabilistic transformation
result.getMetaModel()

[14]:

([0.425629,0.198141] + [0,0.470486] * (-0.29173 + 2.18797 * x0) + [-0.0993209,0.788572] * (0.229006 + 2.32502 * x1) + [-0.246251,0.323295] * (-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) + [-0.192048,-1.84624] * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2) + [0,1.4247] * (-0.207607 - 4.12182 * x1 + 0.977306 * x1^2 + 10.1162 * x1^3) + [-0.476436,0] * ((-0.29173 + 2.18797 * x0) * (0.229006 + 2.32502 * x1)) + [0,-1.02932] * (0.893658 - 1.0145 * x1 - 13.1947 * x1^2 + 2.19409 * x1^3 + 20.4848 * x1^4) + [0,-0.13556] * ((-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) * (0.229006 + 2.32502 * x1)) + [0,-0.281701] * ((-0.29173 + 2.18797 * x0) * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2)) + [-0.0249183,0] * (0.284385 - 7.63723 * x0 - 4.67858 * x0^2 + 78.3356 * x0^3 + 9.59256 * x0^4 - 185.818 * x0^5 - 5.17002 * x0^6 + 120.521 * x0^7) + [0,-0.300647] * (-0.202696 - 7.72688 * x1 + 5.82995 * x1^2 + 87.3021 * x1^3 - 20.8178 * x1^4 - 231.97 * x1^5 + 17.6506 * x1^6 + 168.413 * x1^7) + [0.0755824,0] * ((0.400359 - 4.06203 * x0 - 0.288198 * x0^2 + 7.74638 * x0^3) * (0.229006 + 2.32502 * x1)) + [0,0.340142] * ((-0.29173 + 2.18797 * x0) * (-0.207607 - 4.12182 * x1 + 0.977306 * x1^2 + 10.1162 * x1^3)) + [0.188745,0] * ((-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2)) + [0,0.0346659] * (0.95105 + 2.92377 * x0 - 42.0158 * x0^2 - 17.3514 * x0^3 + 250.399 * x0^4 + 29.114 * x0^5 - 454.203 * x0^6 - 14.8305 * x0^7 + 250.574 * x0^8))o((| y0 = Normal(mu = 0, sigma = 1) -> y0 : Histogram(origin = -1, {w0 = 0.4, h0 = 0.25}, {w1 = 0.4, h1 = 0.25}, {w2 = 0.4, h2 = 0.666667}, {w3 = 0.4, h3 = 1.08333}, {w4 = 0.4, h4 = 0.25})
| y1 = Normal(mu = 0, sigma = 1) -> y1 : Trapezoidal(a = -1, b = -0.530738, c = 0.0480939, d = 1)
)o((class=LinearEvaluation name=Unnamed center=[0,0] constant=[0,0] linear=[[ 1 -0.0108335 ]
[ 0 1.00006 ]])o(| y0 = Histogram(origin = -2.29006, {w0 = 0.805655, h0 = 0.124123}, {w1 = 0.805655, h1 = 0.124123}, {w2 = 0.805655, h2 = 0.330994}, {w3 = 0.805655, h3 = 0.537865}, {w4 = 0.805655, h4 = 0.124123}) -> y0 : Normal(mu = 0, sigma = 1)
| y1 = Trapezoidal(a = -1.75839, b = -0.781359, c = 0.423804, d = 2.40573) -> y1 : Normal(mu = 0, sigma = 1)
)))

[15]:

# Get the projection strategy
algo.getProjectionStrategy()

[15]:

class=ProjectionStrategy implementation=class=LeastSquaresStrategy experiment=class=FixedExperiment name=Unnamed sample=class=Sample name=Normal implementation=class=SampleImplementation name=Normal size=30 dimension=2 description=[X0,X1] data=[[0.608202,-1.26617],[-0.438266,1.20548],[-2.18139,0.350042],[-0.355007,1.43725],[0.810668,0.793156],[-0.470526,0.261018],[-2.29006,-1.28289],[-1.31178,-0.0907838],[0.995793,-0.139453],[-0.560206,0.44549],[0.322925,0.445785],[-1.03808,-0.856712],[0.473617,-0.125498],[0.351418,1.78236],[0.0702074,-0.781366],[-0.721533,-0.241223],[-1.78796,0.40136],[1.36783,1.00434],[0.741548,-0.0436123],[0.539345,0.29995],[0.407717,-0.485112],[-0.382992,-0.752817],[0.257926,1.96876],[-0.671291,1.85579],[0.0521593,0.790446],[0.716353,-0.743622],[0.184356,-1.53073],[0.655027,0.538071],[1.73821,-0.958722],[0.377922,-0.181004]] weights=class=Point name=Unnamed dimension=30 values=[0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333]

OpenTURNS

An Open source initiative for the Treatment of Uncertainties, Risks'N Statistics

Previous topic

Next topic

This Page

Polynomial chaos exploitation¶