Use Kriging as forward model in Bayesian calibration

Python usage
1

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 :

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

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

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

2

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

3

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:

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?

4

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 !

5

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 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])

6

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).

7

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

8

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.

    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 :

    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

9

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:

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)
10

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” ?

11

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.

12

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 :

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 !

13

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:

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:

conditional.getParameterDescription()

Output:

[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).

14

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

1 Like
15

Hi!

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 :slight_smile: .

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

The first example was not clear enough from my point of view, so I tried to improve it in the PR 2368. The proposed example is available at this link. I hope that it will help.

Regards,

Michaël

1 Like
16

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