title: How to specify the prior for scikit-learn's Gaussian process regression?
body: <p>As mentioned <a href="https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.GaussianProcessRegressor.html" rel="noreferrer">here</a>, scikit-learn's Gaussian process regression (GPR) permits "prediction without prior fitting (based on the GP prior)". But I have an idea for what my prior should be (i.e. it should not simply have a mean of zero but perhaps my output, <code>y</code>, scales linearly with my input, <code>X</code>, i.e. <code>y = X</code>). How could I encode this information into GPR?</p>
  
  <p>Below is a working example, but it assumed zero mean for my prior. I <a href="https://blog.dominodatalab.com/fitting-gaussian-process-models-python/" rel="noreferrer">read</a> that "The GaussianProcessRegressor does not allow for the specification of the mean function, always assuming it to be the zero function, highlighting the diminished role of the mean function in calculating the posterior." I believe this is the motivation behind <a href="https://github.com/jmetzen/gp_extras" rel="noreferrer">custom kernels</a> (e.g. heteroscedastic) with variable scales at different <code>X</code>, although I'm still trying to better understand what capability they provide. Are there ways to get around the zero mean prior so that an arbitrary prior can be specified in scikit-learn?</p>
  
  <pre><code>import numpy as np
  from matplotlib import pyplot as plt
  from sklearn.gaussian_process import GaussianProcessRegressor
  from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
  
  def f(x):
      """The function to predict."""
      return 1.5*(1. - np.tanh(100.*(x-0.96))) + 1.5*x*(x-0.95) + 0.4 + 1.5*(1.-x)* np.random.random(x.shape)
  
  # Instantiate a Gaussian Process model
  kernel = C(10.0, (1e-5, 1e5)) * RBF(10.0, (1e-5, 1e5))
  
  X = np.array([0.803,0.827,0.861,0.875,0.892,0.905,
                  0.91,0.92,0.925,0.935,0.941,0.947,0.96,
                  0.974,0.985,0.995,1.0])
  X = np.atleast_2d(X).T
  
  # Observations and noise
  y = f(X).ravel() 
  noise = np.linspace(0.4,0.3,len(X))
  y += noise
  
  # Instantiate a Gaussian Process model
  gp = GaussianProcessRegressor(kernel=kernel, alpha=noise ** 2,
                                n_restarts_optimizer=10)
  # Fit to data using Maximum Likelihood Estimation of the parameters
  gp.fit(X, y)
  
  # Make the prediction on the meshed x-axis (ask for MSE as well)
  x = np.atleast_2d(np.linspace(0.8, 1.02, 1000)).T
  y_pred, sigma = gp.predict(x, return_std=True)
  
  plt.figure() 
  plt.errorbar(X.ravel(), y, noise, fmt='k.', markersize=10, label=u'Observations')
  plt.plot(x, y_pred, 'k-', label=u'Prediction')
  plt.fill(np.concatenate([x, x[::-1]]),
           np.concatenate([y_pred - 1.9600 * sigma,
                          (y_pred + 1.9600 * sigma)[::-1]]),
           alpha=.1, fc='k', ec='None', label='95% confidence interval')
  plt.xlabel('x')
  plt.ylabel('y')
  plt.xlim(0.8, 1.02)
  plt.ylim(0, 5)
  plt.legend(loc='lower left')
  plt.show()
  </code></pre>