Note
Click here to download the full example code
Sample trajectories from a Gaussian Process with correlated outputs¶
A 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 .
import openturns as ot
import openturns.viewer as viewer
from numpy import square
ot.Log.Show(ot.Log.NONE)
Create a Kronecker covariance function¶
First, define the scalar correlation function .
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)
Out:
[[ 1 0.9 ]
[ 0.9 1.5 ]]
Use these ingredients to build the KroneckerCovarianceModel
.
kronecker = ot.KroneckerCovarianceModel(rho, C)
Build a Gaussian Process with Kronecker covariance function¶
Define a 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.setColors(ot.Drawable.BuildDefaultPalette(2))
graph.setLegends(["Y1", "Y2"])
graph.setLegendPosition("topleft")
_ = viewer.View(graph)
![plot kronecker covmodel](../../_images/sphx_glr_plot_kronecker_covmodel_001.png)
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)
_ = viewer.View(graph, square_axes=True)
![Trajectory on the complex plane](../../_images/sphx_glr_plot_kronecker_covmodel_002.png)
Change the correlation between the outputs¶
By setting covariance matrix of the outputs, we have implicitely 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)
Out:
[[ 1 0.9 ]
[ 0.9 1.5 ]]
# Print squared amplitudes
print(square(kronecker.getAmplitude()))
Out:
[1. 1.5]
The diagonal of the correlation matrix is by definition filled with ones.
output_correlation = kronecker.getOutputCorrelation()
print(output_correlation)
Out:
[[ 1 0.734847 ]
[ 0.734847 1 ]]
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)
Out:
[[ 1 0.9 ]
[ 0.9 1 ]]
Changing the output correlation matrix does not change the amplitudes.
kronecker.setOutputCorrelation(output_correlation)
print(square(kronecker.getAmplitude()))
Out:
[1. 1.5]
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.setColors(ot.Drawable.BuildDefaultPalette(2))
graph.setLegends(["Y1", "Y2"])
graph.setLegendPosition("topleft")
_ = viewer.View(graph)
![plot kronecker covmodel](../../_images/sphx_glr_plot_kronecker_covmodel_003.png)
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)
_ = viewer.View(graph, square_axes=True)
![Trajectory on the complex plane](../../_images/sphx_glr_plot_kronecker_covmodel_004.png)
Total running time of the script: ( 0 minutes 0.272 seconds)