# Question about "gp.predict": Gaussian Processes

**URL:** <https://discourse.pymc.io/t/question-about-gp-predict-gaussian-processes/7965>\
**Category:** Questions\
**Created:** [August 29, 2021, 10:29pm UTC](https://discourse.pymc.io/t/question-about-gp-predict-gaussian-processes/7965 "2021-08-29T22:29:21Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![MCruzReyes](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mcruzreyes/32/4392_2.png) [@MCruzReyes](https://discourse.pymc.io/u/MCruzReyes)\
**Post date:** [August 29, 2021, 10:29pm UTC](https://discourse.pymc.io/t/question-about-gp-predict-gaussian-processes/7965/1 "2021-08-29T22:29:21Z")

</div>

I am trying to make predictions using “gp.predict”, however, every time I use “gp.predict” I get a different  
result. The code I am using:

```auto
import matplotlib.pyplot as plt
import pymc3 as pm
import numpy as np
import theano.tensor as tt
import arviz as az
%config InlineBackend.figure_format = 'retina'
az.style.use("arviz-darkgrid")

X = np.array([1,2,3,4])
Y = np.array([1,4,9,16])
Sigma = np.matrix([[1,0.1, 0,0], [0.1, 2,0.2,0.4], [0, 0.2,3,0.4],[0, 0.4,0.4,4]])
le= len(X)

#Reorder the data 
Xg = X.reshape((le,1))

class ConstantMatrix(pm.gp.cov.Covariance):
    def __init__ (self, cov):
        self.cov = tt.as_tensor_variable(cov)
        
    def full(self, X, Xs):
        # covariance matrix known, not explicitly function of X 
        return self.cov

    def diag(self, X):
        return tt.diag(self.cov)

Sigma = ConstantMatrix(np.matrix([[1,0.1, 0,0], [0.1, 2,0.2,0.4], [0, 0.2,3,0.4],[0, 0.4,0.4,4]])) 

with pm.Model() as model:
    ℓ = pm.Gamma("ℓ", alpha= 2, beta=1)
    η = pm.Gamma("η", alpha= 2, beta=1)
    
    cov = η ** 2 * pm.gp.cov.Matern52(1, ℓ)
    gp = pm.gp.Marginal(cov_func=cov)
    
    y_ = gp.marginal_likelihood("y", X=Xg, y=Y, noise=Sigma)

    #MCMC run
    trace = pm.sample(2500, cores=6, tune= 1000,return_inferencedata=True)

# check convergence issues
n_nonconverged = int(np.sum(az.rhat(trace)[["η", "ℓ"]].to_array() > 1.03).values)
print("%i variables MCMC chains appear not to have converged." % n_nonconverged)

```

And the lines of code that give me different results every time I run them are:

```auto
############################################
# Prediction. This part of the code gives different results
X_new = np.linspace(0, np.max(X), 100)[:, None] #Predicted values
mu, var = gp.predict(X_new, point=trace,diag=True)
sd = np.sqrt(var)
############################################

#Plot
fig = plt.figure(figsize=(6, 6))
ax = fig.gca()

#Plot mean and 1σ intervals
plt.plot(X_new, mu, "b")
plt.fill_between(X_new.flatten(), mu + 1 * sd, mu - 1 * sd,facecolor='blue',label= r'Reconstruction 1$\sigma$', alpha=0.4)

#Data
plt.errorbar(X, Y, errors, color='black', fmt='_', label="Data")

#Format
plt.xlim(0, np.max(X))
plt.xlabel('X')
plt.ylabel("Y",fontsize=24)
plt.legend(loc="best", prop={'size': 14})
plt.show()

```

As far as I understand I should get the same result every time I make a prediction, however every time I use the lines of code above I get a different result for the prediction.

I have two questions:

1. What do I need to do to get the same prediction every time I run the code?
2. How can I get the prediction with the values that maximize the “marginal\_likelihood”? I tried to change this line of code

```auto
mu, var = gp.predict(X_new, point=trace,diag=True)

```

by

```auto
mu, var = gp.predict(X_new, point=points,diag=True)

```

where “points” are the values that maximize “marginal\_likelihood”, but it didn’t work.

---

<div class="post-metadata">

**Author:** ![jlindbloom](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jlindbloom/32/5016_2.png) [@jlindbloom](https://discourse.pymc.io/u/jlindbloom)\
**Post date:** [September 2, 2021, 9:08pm UTC](https://discourse.pymc.io/t/question-about-gp-predict-gaussian-processes/7965/2 "2021-09-02T21:08:33Z")

</div>

> 1. What do I need to do to get the same prediction every time I run the code?

Your issue is that you need to pass a single set of model parameters into the kwarg for `point`. I’m not confident if this is actually what’s going on but I think by passing in the entire trace (which contains many samples of the model parameters) it is just randomly selecting a sample to use. For example, if you try

> mu, var = gp.predict(X\_new, point={“η”: 6.0, “ℓ”: 2.1}, diag=True)

you’ll find you get the same results each time. But here I just manually set the values for these two parameters.

> 1. How can I get the prediction with the values that maximize the “marginal\_likelihood”?

You’ll need to run

```Python
with model:
    map_point = pm.find_MAP()

```

in order to get the parameters that maximize the marginal likelihood, and then pass this as your `point` into your code from before.

This is what I get:

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/e/e8ae80047b310c2031583c9e93c01e70c82d5158.png)

---

<div class="post-metadata">

**Author:** ![MCruzReyes](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mcruzreyes/32/4392_2.png) [@MCruzReyes](https://discourse.pymc.io/u/MCruzReyes)\
**Post date:** [September 2, 2021, 11:54pm UTC](https://discourse.pymc.io/t/question-about-gp-predict-gaussian-processes/7965/3 "2021-09-02T23:54:14Z")

</div>

Thank you!