Question about "gp.predict": Gaussian Processes

[MCruzReyes]

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:

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:

############################################
# 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

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

by

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

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

[jlindbloom]

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.

2. How can I get the prediction with the values that maximize the “marginal_likelihood”?

You’ll need to run

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:

[MCruzReyes]

Thank you!