Gaussian process regression: samples from mean model

Questions
1

Hi there!

I have a physical model that I’m fitting against time-series observations. I’ve implemented Gaussian process regression with my (complicated) custom mean model, and I’m hitting a small snag in visualizing the results. I’d really like to draw samples from trace and evaluate the model at those posterior points to generate a noiseless model time-series without the GP regression included. Then, separately, I’d like to visualize the model with the GP.

Ideally I’d like to do something like:

with model: 
    mu_noiseless = model.gp.mean_func(xnew[:, None])

Is this possible with the current PyMC3 API, or do I need to rewrite my mean model with numpy?

Thanks!

2

This might be a bit hacky, but the solution I’m currently working with is to define the Gaussian process as an additive GP where there are two components: a white noise component with my complicated mean model, and a Matern 3/2 component with a zero mean model. This allows me to call the white noise component separate from the Matern component, like so:

mumatern, varmatern = gp_matern.predict(xnew[:, None], point=trace[-1],
                           given={"gp": gp, "X": x[:, None], "y": y, "noise": yerr},
                     diag=True)

muwhite, varwhite = gp_white.predict(xnew[:, None], point=trace[-1],
                           given={"gp": gp, "X": x[:, None], "y": y, "noise": yerr},
                     diag=True)

muwhite then contains the mean model evaluated at trace[-1] and mumatern contains the GP trend. Is this sensible?