# Use Kriging as forward model in Bayesian calibration

**URL:** <https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262>\
**Category:** Python usage\
**Created:** [January 6, 2023, 11:53am UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262 "2023-01-06T11:53:10Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![elieso](https://avatars.discourse-cdn.com/v4/letter/e/f04885/32.png) [@elieso](https://openturns.discourse.group/u/elieso)\
**Post date:** [January 6, 2023, 11:53am UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/1 "2023-01-06T11:53:10Z")

</div>

Hi all,

I am a new user to OpenTURNS but familiar with UQ topics, especially Kriging and Bayesian Calibration through the use of MCMC algorithms.

I am trying to set up a Bayesian calibration problem and use a Kriging model, that i built prior to the calibration step.

I believe that the Kriging model should be transmitted to the calibration algorithm through the use of the function :

```python
mh.setLikelihood(conditional,Yexp,linkFunction,Xdoe)

```

where `linkFunction` is my Kriging object, constructed (in an external script i wrote) with:

```python
krigingMetaModel = result.getMetaModel()

```

It does not work, the `setLikelihood` function outputs an error message about dimensions mismatch.

Maybe i am doing something wrong or i misunderstand the use of this `linkFunction` ?

Thanks by advance for your help !

Best regards,  
Elie SOLAI

---

<div class="post-metadata">

**Author:** ![schueller](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/schueller/32/3_2.png) [@schueller](https://openturns.discourse.group/u/schueller)\
**Post date:** [January 17, 2023, 9:12am UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/2 "2023-01-17T09:12:24Z")

</div>

yes, it can be tricky to get right, the link function must have an input dimension of the chain and the output dimension of the conditional distribution:

[https://openturns.github.io/openturns/latest/user\_manual/\_generated/openturns.RandomVectorMetropolisHastings.html#openturns.RandomVectorMetropolisHastings.setLikelihood](https://openturns.github.io/openturns/latest/user_manual/_generated/openturns.RandomVectorMetropolisHastings.html#openturns.RandomVectorMetropolisHastings.setLikelihood)

---

<div class="post-metadata">

**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\
**Post date:** [January 18, 2023, 9:29am UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/3 "2023-01-18T09:29:14Z")

</div>

Hello Elie, could you perhaps share further details about you are trying to do?

In particular, which function does your Kriging model emulate? I would also like to know more about the other arguments in your code line:

```auto
mh.setLikelihood(conditional,Yexp,linkFunction,Xdoe)

```

- _conditional_ must be a `Distribution` object with parameter(s) \nu. Let f(\cdot | \nu) be its PDF.
- _Yexp_ is the _observations_ argument. Let \boldsymbol y^i denote the i-th observation. Then f(\boldsymbol y^i | \nu) is the PDF at that observation.
- _linkFunction_ (denoted by g) takes the parameter(s) \boldsymbol \theta you want to calibrate as input and outputs a valid value for the \nu parameter(s) of _conditional_. In your case, you say it is a Kriging metamodel, which seems odd to me, because it must take the parameter(s) you want to calibrate as input and output a valid value for \nu, the parameter(s) of _conditional_.
- _Xdoe_ is the _covariates_ argument. It means that _linkFunction_ must be a ParametricFunction allowing _Xdoe[i]_ (denoted below by \boldsymbol x^i) to be its parameter, so that the likelihood is defined as:

L(\boldsymbol \theta) = \prod\_{i=1}^n f(\boldsymbol{y}^i|g\_{\boldsymbol x^i}(\boldsymbol{\theta}))

What I suspect is that your Kriging metamodel really takes \boldsymbol x^i as input (not as parameter) and outputs \boldsymbol y^i, so it is not a valid link function g. Could you share in detail what each argument you provide `setLikelihood` is so I can understand better?

---

<div class="post-metadata">

**Author:** ![elieso](https://avatars.discourse-cdn.com/v4/letter/e/f04885/32.png) [@elieso](https://openturns.discourse.group/u/elieso)\
**Post date:** [January 19, 2023, 10:13am UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/4 "2023-01-19T10:13:14Z")

</div>

Hi,  
many thanks for your detailed answer !  
I’ll try to decompose in details the steps i go through to define the Likelihood function :

Overall in my problem i have 5 uncertain inputs : X = (x\_1,x\_2,...x\_5)  
And one scalar output  
Y = y

My computational model can be denoted \mathcal{M} so i have :  
 Y = \mathcal{M}(X)

My Kriging model denoted \mathcal{M}^K emulates this model \mathcal{M}

First, i take the the Kriging model that i built in an external function (using `krigingMetaModel = result.getMetaModel()` )

> `forwardModel = krigForwardModel`

The inputs of my setLikelihood function are defined in the following :

The conditional distribution is defined with :

> ```
> conditional = ot.Normal()
> 
> ```

Then i define the `linkFunction` using the parametric function. I take :

> `linkFunction = ot.ParametricFunction(forwardModel,[0],[0,0,0,0,0])`

The `[0],[0,0,0,0,0]` are ‘dummy’ values of the output and inputs of my model \mathcal{M}

Then my Yexp is one experimental observation of my QOI:

> Yexp = [[1.05]]

Then the Xdoe is defined with a numpy matrix of size (Nsample,Nparams=5)

I define the random walk MetropolisHastings setting:

> mh\_coll = [  
> ot.RandomWalkMetropolisHastings(priorDistrib, initialState, proposal, [i])  
> for i in range(Nd)  
> ]

Finally the setLikelihood function takes all the parameters I just described above :

> for mh in mh\_coll:  
> mh.setLikelihood(conditional,Yexp,linkFunction,Xdoe)

So from what i am doing and what you explained, I guess I am misunderstanding the role of the linkFunction here ?

Thanks again !

---

<div class="post-metadata">

**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\
**Post date:** [January 21, 2023, 4:44pm UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/5 "2023-01-21T16:44:22Z")

</div>

Thanks for the details! They raise further questions on my part:

> First, i take the the Kriging model that i built in an external function (using `krigingMetaModel = result.getMetaModel()` )
> 
> `forwardModel = krigForwardModel`

I guess you meant `forwardModel = krigingMetaModel`.

> Then i define the `linkFunction` using the parametric function. I take :
> 
> > `linkFunction = ot.ParametricFunction(forwardModel,[0],[0,0,0,0,0])`
> 
> The `[0],[0,0,0,0,0]` are ‘dummy’ values of the output and inputs of my model \mathcal{M}

This is what I don’t understand.The [ParametricFunction — OpenTURNS 1.20 documentation](https://openturns.github.io/openturns/1.20/user_manual/_generated/openturns.ParametricFunction.html) documentation states that the second argument is a list of indices (with type int), not a list of values (with type float): it determines which inputs of `forwardModel` are “frozen” (that is, turned into a parameter). In your case the input of index 0, i.e. the first.

The third argument is the list of values taken by the parameters (e.g., the “frozen” inputs). In your case, since only one input was frozen, the third argument should be a list of length 1, not 5. The line:

`linkFunction = ot.ParametricFunction(forwardModel,[0],[0,0,0,0,0])`

should raise an error, I don’t understand how you were able to run it. I would rather expect something like

`linkFunction = ot.ParametricFunction(forwardModel,[0],[2.3])`

---

<div class="post-metadata">

**Author:** ![elieso](https://avatars.discourse-cdn.com/v4/letter/e/f04885/32.png) [@elieso](https://openturns.discourse.group/u/elieso)\
**Post date:** [January 25, 2023, 4:03pm UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/6 "2023-01-25T16:03:41Z")

</div>

Hi !

Thanks again for these explanations.  
I understand better the use of the ParametricFunction.

This still makes me asking :  
1 / I am wondering if there is an other way to use my kriging model in the linkfunction ? For instance can i use ot.PythonFunction() and give directly the constructed `krigMetaModel` object ?

2/ I still get the error when i give the linkFunction in the mh.setLikelihood :

> mh.setLikelihood(conditional,Yexp,linkFunction)

and still get the dimension mismatch error :

> RuntimeError: InvalidDimensionException : The linkFunction input dimension (0) does not match the dimension of the prior (5).

---

<div class="post-metadata">

**Author:** ![MerlinKeller](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/merlinkeller/32/190_2.png) [@MerlinKeller](https://openturns.discourse.group/u/MerlinKeller)\
**Post date:** [January 26, 2023, 1:15pm UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/7 "2023-01-26T13:15:45Z")

</div>

Hi Elie,

I have been following this thread, and it seems to me the best would be if you could submit us a piece of code that can reproduce the error message, then we could surely identify and solve your bug.

Do you think you could attach such a script to your next post?

Many thanks,  
Merlin

---

<div class="post-metadata">

**Author:** ![elieso](https://avatars.discourse-cdn.com/v4/letter/e/f04885/32.png) [@elieso](https://openturns.discourse.group/u/elieso)\
**Post date:** [January 26, 2023, 5:07pm UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/8 "2023-01-26T17:07:48Z")

</div>

Hi,

Thank you so much for following the thread.

I can attach first the lines where i build the Kriging surrogate model :  
Note that for devpt purposes i just generate here some random data to make the workflow run.

```python
    xdata = np.zeros((Ns,Nd)) # Ns = nb of samples / Nd = nb of input variables
    for i in range(0,Ns): # browse all samples of inputs
        for d in range(0,Nd): # browse nb of variables of the problem 
            xdata[i,d] = i+0.1*d - 0.4**d 

    # outputs / artificial 
    Yqoi = [] 
    for i in range(0,Ns):
        ylocalList = [i*2 + i**2] # artificial data  
        Yqoi.append(ylocalList) 

    # !!! end artificial data generation 

    # KRIGING PARAMETERS 
    print('[main_surrogate] Define Kriging parameters')    
    basis = ot.ConstantBasisFactory(Nd).build()
    covarianceModel = ot.SquaredExponential([0.1] * Nd, [1.0])
    algo = ot.KrigingAlgorithm( xdata , Yqoi , covarianceModel , basis)
    print('[main_surrogate] Run the Kriging construction algorithm') 
    algo.run()
    result = algo.getResult()
    
    # save the constructed Kriging model 
    print('[main_surrogate] Save the Kriging metamodel') 
    krigingMetaModel = result.getMetaModel()

```

Now i can show some parts of the script where i set up the bayesian calibration :

```python
    obsSize = 1 
    Yexp = [[1.05]] # fake experimental value

    # read and get the forward model == the surrogate model 
    forwardModel = krigMetaModel
    # print('[main_bayesian]-DEBUG', forwardModel.getInputDescription()) # debug / 06.01.2023 

    linkFunction = ot.ParametricFunction(forwardModel,[0,1,2,3,4],[2.3,1.2,1.2,1.2,1.4])
    
    # :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    # PRIOR distribution 
    
    marginals = []  
    for i in range(0,Nd): # browse the uncertain variables distribution 
       # set ot ready variable 
        if listUncVar[i].stringDistName == 'Normal':
            locDistrib = ot.Normal(listUncVar[i].vecParams[0], listUncVar[i].vecParams[1])
        elif listUncVar[i].stringDistName == 'Uniform':
            locDistrib = ot.Uniform(listUncVar[i].vecParams[0], listUncVar[i].vecParams[1])
        
        marginals.append( locDistrib ) 
    
    priorDistrib = ot.ComposedDistribution(marginals) # define a multi variate prior based on marginals

    # conditional proposal distribution for the Markov Chain 
    conditional = ot.Normal() # dimension of the output y of the model f: y=f(x)
    proposal = ot.Uniform(-1.0,1.0)
    
    print('[main_bayesian] Compute initial state as mean of distributions') 
    initialState = np.zeros(Nd)
    for i in range(0,Nd):
        if listUncVar[i].stringDistName == 'Normal':
            initialState[i] = listUncVar[i].vecParams[0] # mean of normal distribution      

        elif listUncVar[i].stringDistName == 'Uniform':
            initialState[i] = (listUncVar[i].vecParams[0]+listUncVar[i].vecParams[1])/2.0 # mean of uniform distribution      
    
    mh_coll = [
        ot.RandomWalkMetropolisHastings(priorDistrib,initialState,proposal,[i]) # a random walk for each input  
        for i in range(Nd)
            ]

    for mh in mh_coll:
         mh.setLikelihood(conditional,Yexp,linkFunction)
 
    sampler = ot.Gibbs(mh_coll) 
    sampler.setThinning(1)
    sampler.setBurnIn(2000)
    sampleSize = 1000000
 
    sample = sampler.getSample(sampleSize)
    print('[main_bayesian] Build posterior')
    kernel = ot.KernelSmoothing()
    posterior = kernel.build(sample)
    [mh.getAcceptanceRate() for mh in sampler.getMetropolisHastingsCollection()]

```

Sorry for the big script, but i guess it is necessary to see each step of the bayesian calibration algorithm construction.

Hope this helps to solve the problem !  
Don’t hesitate if you need further lines of code.

Thanks again very much

Elie

---

<div class="post-metadata">

**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\
**Post date:** [January 27, 2023, 10:17am UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/9 "2023-01-27T10:17:04Z")

</div>

Hello, thanks for the script, things are clearer now.

The link function needs to take the Nd-dimensional input and produce the parameters of the distribution of the output. Since that distribution is Normal, it needs to produce its mean and standard deviation.

In the example below (which should work with your script), the mean is the Kriging prediction and the standard deviation is the square root of the Kriging variance:

```python
from math import sqrt

# The link function produces the parameters of the distribution of the observation
def python_link_function(X): # X is the Nd-dimensional input
    """
    The true value of the output is supposed to follow
    a Normal distribution
    centered on the Kriging prediction
    and with variance equal to the Kriging variance.
    """
    mean = forwardModel(X)[0]
    variance = result.getConditionalMarginalVariance(X)
    return [mean, sqrt(variance)] # the parameters of the normal dist are mean and std

linkFunction = ot.PythonFunction(Nd, 2, python_link_function)

```

---

<div class="post-metadata">

**Author:** ![elieso](https://avatars.discourse-cdn.com/v4/letter/e/f04885/32.png) [@elieso](https://openturns.discourse.group/u/elieso)\
**Post date:** [February 2, 2023, 2:15pm UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/10 "2023-02-02T14:15:14Z")

</div>

Ok thanks a lot I get how it works now.

Is it possible to pass the objects such as : `forwardModel` and `results` (obtained after `algo.getResults()`  
) from the Kriging construction script, as inputs of the function “python\_link\_function” ?

---

<div class="post-metadata">

**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\
**Post date:** [February 7, 2023, 10:06am UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/11 "2023-02-07T10:06:07Z")

</div>

Unfortunately no, the “python\_link\_function” is understood by the underlying C++ code as a Function object, which takes a Point as input (or something that can be converted by OpenTURNS to a Point like a list of floats for example) and returns a Point as output. It cannot take other types as input or output.

---

<div class="post-metadata">

**Author:** ![elieso](https://avatars.discourse-cdn.com/v4/letter/e/f04885/32.png) [@elieso](https://openturns.discourse.group/u/elieso)\
**Post date:** [June 14, 2023, 4:21pm UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/12 "2023-06-14T16:21:06Z")

</div>

Hi,

I was slammed by other projects these last few months but i am getting back to this bayesian calibration.

I think i advanced in the understanding and i am close to the right solution with the function you proposed :

```auto
def python_link_function(X): # X is the Nd-dimensional input
    """
    The true value of the output is supposed to follow
    a Normal distribution
    centered on the Kriging prediction
    and with variance equal to the Kriging variance.
    """
    mean = forwardModel(X)[0]
    variance = result.getConditionalMarginalVariance(X)
    return [mean, sqrt(variance)] # the parameters of the normal dist are mean and std

linkFunction = ot.PythonFunction(Nd, 2, python_link_function)

```

However, what if my Kriging model outputs 2 quantities of interest and i want to output from the python\_link\_function, let’s say something like :  
mean\_Y1 = forwardModel(X)[0]  
mean\_Y2 = forwardModel(X)[1]

How can I pass the parameters of the distributions for the 2 outputs of the Kriging (mean,variance) to the linkFunction ?

Thanks again for your (huge) help !

---

<div class="post-metadata">

**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\
**Post date:** [August 16, 2023, 12:01pm UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/13 "2023-08-16T12:01:09Z")

</div>

Hi @elieso and sorry for the delay, I have not had time to browse the forum for a while.

If your Kriging model has output dimension 2, let’s say \boldsymbol{Y}(x) = (Y\_1(x), Y\_2(x))^T where x is the N\_d-dimensional input, then \boldsymbol{Y}(x) follows a bivariate normal distribution.

Therefore your conditional distribution is going to be defined with:

```auto
conditional = ot.Normal(2)

```

And your link function will need to output the parameters of a bivariate normal distribution.  
To see what these parameters are, use the `getParameterDescription` method:

```auto
conditional.getParameterDescription()

```

Output:

```auto
[mean_0,standard_deviation_0,mean_1,standard_deviation_1,R_1_0]

```

which are respectively the mean of Y\_1(x), the standard deviation of Y\_1(x), the mean of Y\_2(x), the standard deviation of Y\_2(x), and the correlation coefficient of Y\_1(x) and Y\_2(x) (with value 0 if they are independent which is by default the case when you do Kriging with multiple outputs).

---

<div class="post-metadata">

**Author:** ![elieso](https://avatars.discourse-cdn.com/v4/letter/e/f04885/32.png) [@elieso](https://openturns.discourse.group/u/elieso)\
**Post date:** [August 23, 2023, 9:55am UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/14 "2023-08-23T09:55:00Z")

</div>

Hi @josephmure !

Thanks very much for your reply.  
It is much more clear to me now, and I just have to say that it works well that way !

Many many thanks !

Elie

---

<div class="post-metadata">

**Author:** ![MichaelBaudin](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/michaelbaudin/32/9_2.png) [@MichaelBaudin](https://openturns.discourse.group/u/MichaelBaudin)\
**Post date:** [August 23, 2023, 4:18pm UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/15 "2023-08-23T16:18:36Z")

</div>

Hi!

> [@elieso](#):
>
> I understand better the use of the ParametricFunction.

Don’t worry: you are not the only one to struggle with these functions. Often, we have issues with calibration features because we do not completely understand this particular class. Actually, it is mandatory to have a clear view on this feature, because everything on top of it only make it more complicated 🙂 .

So let me mention two help pages which better describe how the `ParametricFunction` class works.

- [plot\_parametric\_function](http://openturns.github.io/openturns/latest/auto_functional_modeling/vectorial_functions/plot_parametric_function.html) is an example where we create a function with 2 inputs based on a function with 4 inputs (2 parameters are fixed).
- [plot\_calibration\_quickstart.html#create-the-parametric-function](http://openturns.github.io/openturns/latest/auto_calibration/least_squares_and_gaussian_calibration/plot_calibration_quickstart.html#create-the-parametric-function) is an example where we perform calibration using least squares.

The first example was not clear enough from my point of view, so I tried to improve it in [the PR 2368](https://github.com/openturns/openturns/pull/2368). The proposed example is available [at this link](https://output.circle-artifacts.com/output/job/d5276932-3ff0-4c4f-8c18-ec3c96a17651/artifacts/0/html/auto_functional_modeling/vectorial_functions/plot_parametric_function.html#sphx-glr-auto-functional-modeling-vectorial-functions-plot-parametric-function-py). I hope that it will help.

Regards,

Michaël

---

<div class="post-metadata">

**Author:** ![elieso](https://avatars.discourse-cdn.com/v4/letter/e/f04885/32.png) [@elieso](https://openturns.discourse.group/u/elieso)\
**Post date:** [August 25, 2023, 9:54am UTC](https://openturns.discourse.group/t/use-kriging-as-forward-model-in-bayesian-calibration/262/16 "2023-08-25T09:54:02Z")

</div>

Hi !

Thanks for for the references.  
I love the specific example on the parametric function ‘create a parametric function’.  
It’s crystal clear now !

Thanks !

Best regards,

Elie