'Predict using the Gaussian process regression model
We can also predict based on an unfitted model by using the GP prior.
In addition to the mean of the predictive distribution, also its
standard deviation (return_std=True) or covariance (return_cov=True).
Note that at most one of the two can be requested.
Parameters
X : array-like, shape = (n_samples, n_features)
Query points where the GP is evaluated
return_std : bool, default: False
If True, the standard-deviation of the predictive distribution at
the query points is returned along with the mean.
return_cov : bool, default: False
If True, the covariance of the joint predictive distribution at
the query points is returned along with the mean
Returns
y_mean : array, shape = (n_samples, [n_output_dims])
Mean of predictive distribution a query points
y_std : array, shape = (n_samples,), optional
Standard deviation of predictive distribution at query points.
Only returned when return_std is True.
y_cov : array, shape = (n_samples, n_samples), optional
Covariance of joint predictive distribution a query points.
Only returned when return_cov is True.'

def predict(self, X, return_std=False, return_cov=False):

    if (return_std and return_cov):
        raise RuntimeError('Not    returning    standard    deviation    of    predictions    when    returning    full    covariance.')
    X = check_array(X)
    if (not hasattr(self, 'X_train_')):
        if (self.kernel is None):
            kernel = (C(1.0, constant_value_bounds='fixed') * RBF(1.0, length_scale_bounds='fixed'))
        else:
            kernel = self.kernel
        y_mean = np.zeros(X.shape[0])
        if return_cov:
            y_cov = kernel(X)
            return (y_mean, y_cov)
        elif return_std:
            y_var = kernel.diag(X)
            return (y_mean, np.sqrt(y_var))
        else:
            return y_mean
    else:
        K_trans = self.kernel_(X, self.X_train_)
        y_mean = K_trans.dot(self.alpha_)
        y_mean = (self._y_train_mean + y_mean)
        if return_cov:
            v = cho_solve((self.L_, True), K_trans.T)
            y_cov = (self.kernel_(X) - K_trans.dot(v))
            return (y_mean, y_cov)
        elif return_std:
            L_inv = solve_triangular(self.L_.T, np.eye(self.L_.shape[0]))
            K_inv = L_inv.dot(L_inv.T)
            y_var = self.kernel_.diag(X)
            y_var -= np.einsum('ij,ij->i', np.dot(K_trans, K_inv), K_trans)
            y_var_negative = (y_var < 0)
            if np.any(y_var_negative):
                warnings.warn('Predicted    variances    smaller    than    0.    Setting    those    variances    to    0.')
                y_var[y_var_negative] = 0.0
            return (y_mean, np.sqrt(y_var))
        else:
            return y_mean
