KrigingAlgorithm¶

(Source code, png, hires.png, pdf)

class KrigingAlgorithm(*args)¶

Kriging algorithm.

Refer to Kriging.

Available constructors:

KrigingAlgorithm(inputSample, outputSample, covarianceModel, basis)

KrigingAlgorithm(inputSample, outputSample, covarianceModel, basisCollection)

Parameters

inputSample, outputSample2-d sequence of float

The samples $(\vect{x}_k)_{1 \leq k \leq N} \in \Rset^d$ and $(\vect{y}_k)_{1 \leq k \leq N}\in \Rset^p$ upon which the meta-model is built.

covarianceModelCovarianceModel

Covariance model used for the underlying Gaussian process assumption.

basisBasis

Functional basis to estimate the trend (universal kriging): $(\varphi_j)_{1 \leq j \leq n_1}: \Rset^d \rightarrow \Rset$ .

If $p>1$ , the same basis is used for each marginal output.

basisCollectionsequence of Basis

Collection of $p$ functional basis: one basis for each marginal output: $\left[(\varphi_j^1)_{1 \leq j \leq n_1}, \dots, (\varphi_j^p)_{1 \leq j \leq n_p}\right]$ . If the sequence is empty, no trend coefficient is estimated (simple kriging).

Notes

We suppose we have a sample $(\vect{x}_k, \vect{y}_k)_{1 \leq k \leq N}$ where $\vect{y}_k = \cM(\vect{x}_k)$ for all k, with $\cM:\Rset^d \mapsto \Rset^p$ the model.

The meta model Kriging is based on the same principles as those of the general linear model: it assumes that the sample $(\vect{y}_k)_{1 \leq k \leq N}$ is considered as the trace of a Gaussian process $\vect{Y}(\omega, \vect{x})$ on $(\vect{x}_k)_{1 \leq k \leq N}$ . The Gaussian process $\vect{Y}(\omega, \vect{x})$ is defined by:

(1)¶ $\vect{Y}(\omega, \vect{x}) = \vect{\mu}(\vect{x}) + W(\omega, \vect{x})$

where:

$\vect{\mu}(\vect{x}) = \left( \begin{array}{l} \mu_1(\vect{x}) \\ \dots \\ \mu_p(\vect{x}) \end{array} \right)$

with $\mu_l(\vect{x}) = \sum_{j=1}^{n_l} \beta_j^l \varphi_j^l(\vect{x})$ and $\varphi_j^l: \Rset^d \rightarrow \Rset$ the trend functions.

$W$ is a Gaussian process of dimension p with zero mean and covariance function $C = C(\vect{\theta}, \vect{\sigma}, \mat{R}, \vect{\lambda})$ (see CovarianceModel for the notations).

The estimation of the all parameters (the trend coefficients $\beta_j^l$ , the scale $\vect{\theta}$ and the amplitude $\vect{\sigma}$ ) are made by the GeneralLinearModelAlgorithm class.

The Kriging algorithm makes the general linear model interpolating on the input samples. The Kriging meta model $\tilde{\cM}$ is defined by:

$\tilde{\cM}(\vect{x}) = \vect{\mu}(\vect{x}) + \Expect{\vect{Y}(\omega, \vect{x})\, | \, \cC}$

where $\cC$ is the condition $\vect{Y}(\omega, \vect{x}_k) = \vect{y}_k$ for each $k \in [1, N]$ .

(1) writes:

$\tilde{\cM}(\vect{x}) = \vect{\mu}(\vect{x}) + \Cov{\vect{Y}(\omega, \vect{x}), (\vect{Y}(\omega, \vect{x}_1), \dots, \vect{Y}(\omega, \vect{x}_N))} \vect{\gamma}$

where $\Cov{\vect{Y}(\omega, \vect{x}), (\vect{Y}(\omega, \vect{x}_1), \dots, \vect{Y}(\omega, \vect{x}_N))} = \left( \mat{C}( \vect{x}, \vect{x}_1) | \dots | \mat{C}( \vect{x}, \vect{x}_N) \right)$ is a matrix in $\cM_{p,NP}(\Rset)$ and $\vect{\gamma} = \mat{C}^{-1}(\vect{y}-\vect{m})$ .

A known centered gaussian observation noise $\epsilon_k$ can be taken into account with setNoise():

$\hat{\vect{y}}_k = \vect{y}_k + \epsilon_k, \epsilon_k \sim \mathcal{N}(0, \tau_k^2)$

Examples

Create the model $\cM: \Rset \mapsto \Rset$ and the samples:

>>> import openturns as ot
>>> f = ot.SymbolicFunction(['x'],  ['x * sin(x)'])
>>> sampleX = [[1.0], [2.0], [3.0], [4.0], [5.0], [6.0], [7.0], [8.0]]
>>> sampleY = f(sampleX)

Create the algorithm:

>>> basis = ot.Basis([ot.SymbolicFunction(['x'], ['x']), ot.SymbolicFunction(['x'], ['x^2'])])
>>> covarianceModel = ot.SquaredExponential([1.0])
>>> covarianceModel.setActiveParameter([])
>>> algo = ot.KrigingAlgorithm(sampleX, sampleY, covarianceModel, basis)
>>> algo.run()

Get the resulting meta model:

>>> result = algo.getResult()
>>> metamodel = result.getMetaModel()

Methods

`BuildDistribution`(inputSample)	Recover the distribution, with metamodel performance in mind.
`getClassName`()	Accessor to the object’s name.
`getDistribution`()	Accessor to the joint probability density function of the physical input vector.
`getId`()	Accessor to the object’s id.
`getInputSample`()	Accessor to the input sample.
`getMethod`()	Linear algebra method accessor.
`getName`()	Accessor to the object’s name.
`getNoise`()	Observation noise variance accessor.
`getOptimizationAlgorithm`()	Accessor to solver used to optimize the covariance model parameters.
`getOptimizationBounds`()	Accessor to the optimization bounds.
`getOptimizeParameters`()	Accessor to the covariance model parameters optimization flag.
`getOutputSample`()	Accessor to the output sample.
`getReducedLogLikelihoodFunction`()	Accessor to the reduced log-likelihood function that writes as argument of the covariance’s model parameters.
`getResult`()	Get the results of the metamodel computation.
`getShadowedId`()	Accessor to the object’s shadowed id.
`getVisibility`()	Accessor to the object’s visibility state.
`hasName`()	Test if the object is named.
`hasVisibleName`()	Test if the object has a distinguishable name.
`run`()	Compute the response surface.
`setDistribution`(distribution)	Accessor to the joint probability density function of the physical input vector.
`setMethod`(method)	Linear algebra method set accessor.
`setName`(name)	Accessor to the object’s name.
`setNoise`(noise)	Observation noise variance accessor.
`setOptimizationAlgorithm`(solver)	Accessor to the solver used to optimize the covariance model parameters.
`setOptimizationBounds`(optimizationBounds)	Accessor to the optimization bounds.
`setOptimizeParameters`(optimizeParameters)	Accessor to the covariance model parameters optimization flag.
`setShadowedId`(id)	Accessor to the object’s shadowed id.
`setVisibility`(visible)	Accessor to the object’s visibility state.

__init__(*args)¶: Initialize self. See help(type(self)) for accurate signature.

static BuildDistribution(inputSample)¶

Recover the distribution, with metamodel performance in mind.

For each marginal, find the best 1-d continuous parametric model else fallback to the use of a nonparametric one.

The selection is done as follow:

We start with a list of all parametric models (all factories)

For each model, we estimate its parameters if feasible.

We check then if model is valid, ie if its Kolmogorov score exceeds a threshold fixed in the MetaModelAlgorithm-PValueThreshold ResourceMap key. Default value is 5%

We sort all valid models and return the one with the optimal criterion.

For the last step, the criterion might be BIC, AIC or AICC. The specification of the criterion is done through the MetaModelAlgorithm-ModelSelectionCriterion ResourceMap key. Default value is fixed to BIC. Note that if there is no valid candidate, we estimate a non-parametric model (KernelSmoothing or Histogram). The MetaModelAlgorithm-NonParametricModel ResourceMap key allows selecting the preferred one. Default value is Histogram

One each marginal is estimated, we use the Spearman independence test on each component pair to decide whether an independent copula. In case of non independence, we rely on a NormalCopula.

Parameters

sampleSample: Input sample.

Returns

distributionDistribution: Input distribution.

getClassName()¶

Accessor to the object’s name.

Returns

class_namestr: The object class name (object.__class__.__name__).

getDistribution()¶

Accessor to the joint probability density function of the physical input vector.

Returns

distributionDistribution: Joint probability density function of the physical input vector.

getId()¶

Accessor to the object’s id.

Returns

idint: Internal unique identifier.

getInputSample()¶

Accessor to the input sample.

Returns

inputSampleSample: The input sample $(\vect{x}_k)_{1 \leq k \leq N}$ .

getMethod()¶

Linear algebra method accessor.

Returns

methodstr: Used linear algebra method.

getName()¶

Accessor to the object’s name.

Returns

namestr: The name of the object.

getNoise()¶

Observation noise variance accessor.

Returns

noisesequence of positive float: The noise variance $\tau_k^2$ of each output value.

getOptimizationAlgorithm()¶

Accessor to solver used to optimize the covariance model parameters.

Returns

algorithmOptimizationAlgorithm: Solver used to optimize the covariance model parameters.

getOptimizationBounds()¶

Accessor to the optimization bounds.

Returns

problemInterval: The bounds used for numerical optimization of the likelihood.

getOptimizeParameters()¶

Accessor to the covariance model parameters optimization flag.

Returns

optimizeParametersbool: Whether to optimize the covariance model parameters.

getOutputSample()¶

Accessor to the output sample.

Returns

outputSampleSample: The output sample $(\vect{y}_k)_{1 \leq k \leq N}$ .

getReducedLogLikelihoodFunction()¶

Accessor to the reduced log-likelihood function that writes as argument of the covariance’s model parameters.

Returns

reducedLogLikelihoodFunction: The potentially reduced log-likelihood function.

Notes

We use the same notations as in CovarianceModel and GeneralLinearModelAlgorithm : $\vect{\theta}$ refers to the scale parameters and $\vect{\sigma}$ the amplitude. We can consider three situtations here:

Output dimension is $\geq 2$ . In that case, we get the full log-likelihood function $\mathcal{L}(\vect{\theta}, \vect{\sigma})$ .

Output dimension is 1 and the GeneralLinearModelAlgorithm-UseAnalyticalAmplitudeEstimate key of ResourceMap is set to True. The amplitude parameter of the covariance model $\vect{\theta}$ is in the active set of parameters and thus we get the reduced log-likelihood function $\mathcal{L}(\vect{\theta})$ .

Output dimension is 1 and the GeneralLinearModelAlgorithm-UseAnalyticalAmplitudeEstimate key of ResourceMap is set to False. In that case, we get the full log-likelihood $\mathcal{L}(\vect{\theta}, \vect{\sigma})$ .

The reduced log-likelihood function may be useful for some pre/postprocessing: vizualisation of the maximizer, use of an external optimizers to maximize the reduced log-likelihood etc.

Examples

Create the model $\cM: \Rset \mapsto \Rset$ and the samples:

>>> import openturns as ot
>>> f = ot.SymbolicFunction(['x0'], ['x0 * sin(x0)'])
>>> inputSample = ot.Sample([[1.0], [3.0], [5.0], [6.0], [7.0], [8.0]])
>>> outputSample = f(inputSample)

Create the algorithm:

>>> basis = ot.ConstantBasisFactory().build()
>>> covarianceModel = ot.SquaredExponential(1)
>>> algo = ot.KrigingAlgorithm(inputSample, outputSample, covarianceModel, basis)
>>> algo.run()

Get the reduced log-likelihood function:

>>> reducedLogLikelihoodFunction = algo.getReducedLogLikelihoodFunction()

getResult()¶

Get the results of the metamodel computation.

Returns

resultKrigingResult: Structure containing all the results obtained after computation and created by the method run().

getShadowedId()¶

Accessor to the object’s shadowed id.

Returns

idint: Internal unique identifier.

getVisibility()¶

Accessor to the object’s visibility state.

Returns

visiblebool: Visibility flag.

hasName()¶

Test if the object is named.

Returns

hasNamebool: True if the name is not empty.

hasVisibleName()¶

Test if the object has a distinguishable name.

Returns

hasVisibleNamebool: True if the name is not empty and not the default one.

run()¶

Compute the response surface.

Notes

It computes the kriging response surface and creates a KrigingResult structure containing all the results.

setDistribution(distribution)¶

Accessor to the joint probability density function of the physical input vector.

Parameters

distributionDistribution: Joint probability density function of the physical input vector.

setMethod(method)¶

Linear algebra method set accessor.

Parameters

methodstr: Used linear algebra method. Value should be LAPACK or HMAT

Notes

The setter update the implementation and require new evaluation. We might also use the ResourceMap key to set the method when instantiating the algorithm. For that purpose, we can use ResourceMap.SetAsString(GeneralLinearModelAlgorithm-LinearAlgebra, key) with key being HMAT or LAPACK.

setName(name)¶

Accessor to the object’s name.

Parameters

namestr: The name of the object.

setNoise(noise)¶

Observation noise variance accessor.

Parameters

noisesequence of positive float: The noise variance $\tau_k^2$ of each output value.

setOptimizationAlgorithm(solver)¶

Accessor to the solver used to optimize the covariance model parameters.

Parameters

algorithmOptimizationAlgorithm: Solver used to optimize the covariance model parameters.

Examples

Create the model $\cM: \Rset \mapsto \Rset$ and the samples:

>>> import openturns as ot
>>> input_data = ot.Uniform(-1.0, 2.0).getSample(10)
>>> model = ot.SymbolicFunction(['x'], ['x-1+sin(pi_*x/(1+0.25*x^2))'])
>>> output_data = model(input_data)

Create the Kriging algorithm with the optimizer option:

>>> basis = ot.Basis([ot.SymbolicFunction(['x'], ['0.0'])])
>>> thetaInit = 1.0
>>> covariance = ot.GeneralizedExponential([thetaInit], 2.0)
>>> bounds = ot.Interval(1e-2,1e2)
>>> algo = ot.KrigingAlgorithm(input_data, output_data, covariance, basis)
>>> algo.setOptimizationBounds(bounds)

setOptimizationBounds(optimizationBounds)¶

Accessor to the optimization bounds.

Parameters

boundsInterval: The bounds used for numerical optimization of the likelihood.

Notes

See GeneralLinearModelAlgorithm class for more details, particularly setOptimizationBounds().

setOptimizeParameters(optimizeParameters)¶

Accessor to the covariance model parameters optimization flag.

Parameters

optimizeParametersbool: Whether to optimize the covariance model parameters.

setShadowedId(id)¶

Accessor to the object’s shadowed id.

Parameters

idint: Internal unique identifier.

setVisibility(visible)¶

Accessor to the object’s visibility state.

Parameters

visiblebool: Visibility flag.

OpenTURNS

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

Previous topic

Next topic

This Page

KrigingAlgorithm¶