@trot Sorry it took me a bit to get back to you.
I think your use case is a bit different than mine but we can try talking through it. So, from what you wrote, you have fit a GP model in GPy and would like to sample using that model in some way in pymc3 either by using it to represent some quantity of interest or to model some complicated likelihood function directly is that correct?
If you are trying to use a GPy fit directly in a pymc3 model I think you will run into some difficulties because pymc3 requires the gradients w.r.t. the parameters in order to use NUTS which I don’t believe GPy can give you. If you still want to do that you will need to wrap the function using something like this approach for a black-box likelihood function.
Note: If you use this approach you will have to use sub-optimal samplers!
If you are trying to re-create GPy in pymc3 ( which I have done to some extent ) you can try and implement the basic Gaussian process model directly using Theano types. If that’s the case, here is the textbook I used when I was doing this myself.
As an alternative ( and what I would suggest you consider ):
You could use the python Surrogate Modeling Toolbox which has all kinds of surrogate models and, most crucially, also includes their gradients. You can wrap these functions using the black-box tutorial, and a bit of finagling, and have a collection of surrogate models that work directly in pymc3.
Hopefully that is helpful!