Hello,
I’m also very new here, but I can direct you to my post : Bayesian calibration on GP model.
In hope that it yields some useful information 
In short, you may want to build your GP within PyMC, using marginal_likelihood for specifying x and y values. I don’t know if sklearn GP can be incorporated as is in PyMC ecosystem (notably for gradient computation), but you can perhaps relapse to black-box likelihoods (Using a “black box” likelihood function — PyMC3 3.11.5 documentation)