# Gaussian process regression: samples from mean model

**URL:** <https://discourse.pymc.io/t/gaussian-process-regression-samples-from-mean-model/5695>\
**Category:** Questions\
**Created:** [August 25, 2020, 8:35am UTC](https://discourse.pymc.io/t/gaussian-process-regression-samples-from-mean-model/5695 "2020-08-25T08:35:14Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![bmmorris](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bmmorris/32/3851_2.png) [@bmmorris](https://discourse.pymc.io/u/bmmorris)\
**Post date:** [August 25, 2020, 8:35am UTC](https://discourse.pymc.io/t/gaussian-process-regression-samples-from-mean-model/5695/1 "2020-08-25T08:35:14Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![bmmorris](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bmmorris/32/3851_2.png) [@bmmorris](https://discourse.pymc.io/u/bmmorris)\
**Post date:** [August 25, 2020, 9:10am UTC](https://discourse.pymc.io/t/gaussian-process-regression-samples-from-mean-model/5695/2 "2020-08-25T09:10:50Z")

</div>

This might be a bit hacky, but the solution I’m currently working with is to define the Gaussian process as an [additive GP](https://docs.pymc.io/Gaussian_Processes.html#additive-gps) 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?