Bayesian calibration on GP model

Questions v3
1

I just started discovering pyMC few days ago, so please forgive me if my question is really basic, I couldn’t find or understand the answer yet.

My objective is very straightforward, I have X and y values from a given code, on which I can build a GP model. I usually use GPy, and get maximum likelihood hyperparameters, yet why not remain fully Bayesian in pyMC. Once this GP is built (fixed hyperparams and conditioned), I wish to solve an inverse problem in a Bayesian manner. I have an experimental y_exp value, and would like to get a posterior distribution of x_exp. In short, with my current GP model accuracy, how certain am I on the calibrated x_exp.

This is were I struggle, Can this be done seemlessly in pyMC3, with an “observed” argument somewhere (in gp.conditional or predict ?) ? Can this be solved with gradient-based samplers ? Or should I fall back to a black-box likelihood that I would calculate myself ?

Thanks a lot in advance

2

From your description I think it should be possible. This thread may be helpful to get you started, Bayesian model calibration with Gaussian Process - #18 by kejzlarv

1 Like
3

Thank you for the link. I had already seen this post but couldn’t quite find what was missing in my understanding. After going over it again, I think what I was looking for was the gp.marginal_likelihood function, that “tells” that some values were observed. Giving a free variable to the X argument would then allow to recover it, which is my objective.

However, I now struggle on shapes and theano conversion. In my toy example, I have a 1D GP, with a scalar observed value and a scalar parameter theta to recover. It seems that I can’t find how to properly give theta to the gp.marginal_likelihood function. Either the number of dimensions is too low or theano doesn’t want object typed variables.

Here is the example:

import numpy as np
import pymc3 as pm
import theano.tensor as tt

def f(x):
    return np.sin(10*x)

theta_train = np.linspace(0,1,7)
theta_true = 0.8

f_train = f(theta_train)
f_true = f(theta_true)

with pm.Model() as model:
    
    theta = pm.Uniform("theta", lower = 0, upper = 1)
    sigma = pm.HalfNormal("sigma", sd = 1)
    l_scale = pm.Gamma("l_scale", alpha= 1, beta = 5, shape = 2)
    s_n = pm.HalfNormal("s_n", sd = 1)
    
    cov = pm.gp.cov.Matern52(1, l_scale)
    
    gp = pm.gp.Marginal(cov_func = cov)
    
    theta_all = np.concatenate((theta_train,[theta]))
    f_all = np.concatenate((f_train,[f_true]))
    
    y1 = gp.marginal_likelihood("y1",
                                X=theta_all.reshape(-1,1),
                                y=f_all.reshape(-1,1),
                                noise=s_n)

Again I’m sorry if the answer only relies on basic theano shape handling.
Thank you for your lights !

4

Hi, if you new to PyMC, it may be better to use PyMC v4 with Aesara. The new PyMC version can handle shape better, and there are a lot of improvement in Aesara compared to Theano.

For the shape handling in PyMC3, you can check this notebook: gp_regression/spawning_salmon.ipynb at master · fonnesbeck/gp_regression · GitHub

I converted the above notebook into PyMC v4 here: gp_experiments/gp_salmon_pymc_v4.ipynb at main · danhphan/gp_experiments · GitHub

For your exmaple, in this line of codes:

    y1 = gp.marginal_likelihood("y1",
                                X=theta_all.reshape(-1,1),
                                y=f_all.reshape(-1,1),
                                noise=s_n)

Maybe X should be shape (n,1) while y should be shape (n,) but not (n,1).

2 Likes
5

Thank you for your answer, I will gladly have a look at these notebooks. I will try PyMC v4 when possible, to compare with PyMC3.

Yet in my case, I mostly face this error:

TypeError: Unsupported dtype for TensorType: object

which, in my understanding, is caused by the fact that I put a pymc3 free variable in the X argument.

6

Unfortunately, I haven’t been able to solve my issue yet. Has anyone already given pyMC variables to X (concatenated with known values) in the marginal_likelihood method in pyMC3 ?
Many thanks !

7

Hi there MickaRiv,

I was able to get it to run with a few minor tweaks to your code. I only have PyMC v4 installed, so I switched to aesara and there were a few small syntax differences. Apart from that, I fiddled with the reshapes a bit and used pm.math.concatenate for the theta concatenation. I think the latter was the key thing needed to get it to work.

import numpy as np
import pymc as pm
import aesara.tensor as at

def f(x):
    return np.sin(10*x)

theta_train = np.linspace(0,1,7)
theta_true = 0.8

f_train = f(theta_train)
f_true = f(theta_true)

with pm.Model() as model:
    
    theta = pm.Uniform("theta", lower = 0, upper = 1)
    sigma = pm.HalfNormal("sigma", 1)
    l_scale = pm.Gamma("l_scale", alpha= 1, beta = 5, shape = 2)
    s_n = pm.HalfNormal("s_n", 1)
    
    cov = pm.gp.cov.Matern52(1, l_scale)
    
    gp = pm.gp.Marginal(cov_func = cov)
    
    theta_all = pm.math.concatenate((theta_train,[theta]))
    f_all = np.concatenate((f_train,[f_true]))
    
    y1 = gp.marginal_likelihood("y1",
                                X=theta_all.reshape((-1,1)),
                                y=f_all.reshape((-1,)),
                                noise=s_n)

Hope this works for you and is helpful!

3 Likes
8

Hi Martin_Ingram,
Thanks a lot for your answer. I took some time to test this on my side, but it indeed works like a charm.
pm.math.concatenate does seem to be the key here.
It also works properly with pyMC3.

Thanks again !

2 Likes