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
|
"""
Sample trajectories from a Gaussian Process with correlated outputs
===================================================================
A :class:`~openturns.KroneckerCovarianceModel` takes a covariance function
with 1-d output and makes its output multidimensional.
In this example, we use one to define a Gaussian Process with 2 correlated scalar outputs.
For that purpose, a covariance matrix for the outputs is needed
in addition to a scalar correlation function :math:`\\rho`.
"""
# %%
import openturns as ot
import openturns.viewer as otv
from numpy import square
# %%
# Create a Kronecker covariance function
# --------------------------------------
#
# First, define the scalar correlation function :math:`\rho`.
# %%
theta = [2.0]
rho = ot.MaternModel(theta, 1.5)
# %%
# Second, define the covariance matrix of the outputs.
# %%
C = ot.CovarianceMatrix(2)
C[0, 0] = 1.0
C[1, 1] = 1.5
C[1, 0] = 0.9
print(C)
# %%
# Use these ingredients to build the :class:`~openturns.KroneckerCovarianceModel`.
# %%
kronecker = ot.KroneckerCovarianceModel(rho, C)
# %%
# Build a Gaussian Process with Kronecker covariance function
# -----------------------------------------------------------
#
# Define a :class:`~openturns.GaussianProcess` with null trend using this covariance function.
# %%
gp = ot.GaussianProcess(kronecker, ot.RegularGrid(0.0, 0.1, 100))
# %%
# Sample and draw a realization of the Gaussian process.
# %%
ot.RandomGenerator.SetSeed(5)
realization = gp.getRealization()
graph = realization.drawMarginal(0)
graph.add(realization.drawMarginal(1))
graph.setYTitle("")
graph.setTitle("")
graph.setLegends(["Y1", "Y2"])
graph.setLegendPosition("upper left")
_ = otv.View(graph)
# %%
# Draw the trajectory on the complex plane.
# %%
graph = realization.draw()
graph.setXTitle("Real part")
graph.setYTitle("Imaginary part")
graph.setTitle("Trajectory on the complex plane")
diagonal = ot.Curve([-1.5, 1.5], [-1.5, 1.5])
diagonal.setLineStyle("dotdash")
diagonal.setColor("grey")
graph.add(diagonal)
_ = otv.View(graph, square_axes=True)
# %%
# Change the correlation between the outputs
# ------------------------------------------
#
# By setting covariance matrix of the outputs, we have implicitly set the
# amplitudes and the correlation matrix of the Kronecker covariance function.
#
# The amplitudes are the square roots of the diagonal elements
# of the covariance matrix.
# %%
# Recall C
print(C)
# %%
# Print squared amplitudes
print(square(kronecker.getAmplitude()))
# %%
# The diagonal of the correlation matrix is by definition filled with ones.
# %%
output_correlation = kronecker.getOutputCorrelation()
print(output_correlation)
# %%
# Since the correlation matrix is symmetric
# and its diagonal necessarily contains ones,
# we only need to change its upper right (or bottom left) element.
# %%
output_correlation[0, 1] = 0.9
print(output_correlation)
# %%
# Changing the output correlation matrix does not change the amplitudes.
# %%
kronecker.setOutputCorrelation(output_correlation)
print(square(kronecker.getAmplitude()))
# %%
# Let us resample a trajectory.
#
# To show the effect ot the output correlation change,
# we use the same random seed in order to get a comparable trajectory.
# %%
gp = ot.GaussianProcess(kronecker, ot.RegularGrid(0.0, 0.1, 100))
ot.RandomGenerator.SetSeed(5)
realization = gp.getRealization()
graph = realization.drawMarginal(0)
graph.add(realization.drawMarginal(1))
graph.setYTitle("")
graph.setTitle("")
graph.setLegends(["Y1", "Y2"])
graph.setLegendPosition("upper left")
_ = otv.View(graph)
# %%
graph = realization.draw()
graph.setXTitle("Real part")
graph.setYTitle("Imaginary part")
graph.setTitle("Trajectory on the complex plane")
diagonal = ot.Curve([-1.5, 1.5], [-1.5, 1.5])
diagonal.setLineStyle("dotdash")
diagonal.setColor("grey")
graph.add(diagonal)
_ = otv.View(graph, square_axes=True)
# %%
# Display all figures
otv.View.ShowAll()
|
