Interpolating an expensive physical model on a pre-computed grid

I’m having trouble finding the analog in pymc3 for interpolating a several-dimensional grid. Please excuse if I abuse notation or terminology, I’ve tried to boil this question down to its most basic…

I have some observables a, b, c, etc., which are determined by some unobserved quantities p, q, and r, something like:

a = F(p,q,r) + noise
b = G(p,q,r) + noise
c = H(p,q,r) + noise
...

and in this case, I’d like to infer p, q, and r.

The functions F, G, and H are in principle deterministic and relatively smooth, but very expensive to compute, so it’s not practical to include them explicitly in a model and sample away. Instead, it is standard practice in my field to pre-compute those functions on regular grids of their arguments (p, q, and r), and interpolate between the points to make predictions of a, b, and c. Ordinarily, I’d interpolate in numpy, but given the rigamarole of wrapping with theano.as_op, I wonder if there’s a different/simpler way.

Specifically, what about thinking of F, G, and H as Gaussian Processes? My intuition is that under the assumption that my pre-computation on the grid captures most of the meaningful variation in the physical model, this approach could work like an interpolation, albeit an explicitly differentiable one, which means I can use HMC. Has anyone done this, or does anyone have thoughts?