def predict(self, Xnew, full_cov=False, Y_metadata=None, kern=None,
                likelihood=None, include_likelihood=True):
        """
        Predict the function(s) at the new point(s) Xnew. This includes the
        likelihood variance added to the predicted underlying function
        (usually referred to as f).

        In order to predict without adding in the likelihood give
        `include_likelihood=False`, or refer to self.predict_noiseless().

        :param Xnew: The points at which to make a prediction
        :type Xnew: np.ndarray (Nnew x self.input_dim)
        :param full_cov: whether to return the full covariance matrix, or just
                         the diagonal
        :type full_cov: bool
        :param Y_metadata: metadata about the predicting point to pass to the
                           likelihood
        :param kern: The kernel to use for prediction (defaults to the model
                     kern). this is useful for examining e.g. subprocesses.
        :param include_likelihood: Whether or not to add likelihood noise to
                                   the predicted underlying latent function f.
        :type include_likelihood: bool

        :returns: (mean, var):
            mean: posterior mean, a Numpy array, Nnew x self.input_dim
            var: posterior variance, a Numpy array, Nnew x 1 if full_cov=False,
                 Nnew x Nnew otherwise

            If full_cov and self.input_dim > 1, the return shape of var is
            Nnew x Nnew x self.input_dim. If self.input_dim == 1, the return
            shape is Nnew x Nnew. This is to allow for different normalizations
            of the output dimensions.

        Note: If you want the predictive quantiles (e.g. 95% confidence
        interval) use :py:func:"~GPy.core.gp.GP.predict_quantiles".
        """

        # Predict the latent function values
        mean, var = self._raw_predict(Xnew, full_cov=full_cov, kern=kern)

        if include_likelihood:
            # now push through likelihood
            if likelihood is None:
                likelihood = self.likelihood
            mean, var = likelihood.predictive_values(mean, var, full_cov,
                                                     Y_metadata=Y_metadata)

        if self.normalizer is not None:
            mean = self.normalizer.inverse_mean(mean)

            # We need to create 3d array for the full covariance matrix with
            # multiple outputs.
            if full_cov & (mean.shape[1] > 1):
                var = self.normalizer.inverse_covariance(var)
            else:
                var = self.normalizer.inverse_variance(var)

        return mean, var
Predict the function(s) at the new point(s) Xnew. This includes the
        likelihood variance added to the predicted underlying function
        (usually referred to as f).

        In order to predict without adding in the likelihood give
        `include_likelihood=False`, or refer to self.predict_noiseless().

        :param Xnew: The points at which to make a prediction
        :type Xnew: np.ndarray (Nnew x self.input_dim)
        :param full_cov: whether to return the full covariance matrix, or just
                         the diagonal
        :type full_cov: bool
        :param Y_metadata: metadata about the predicting point to pass to the
                           likelihood
        :param kern: The kernel to use for prediction (defaults to the model
                     kern). this is useful for examining e.g. subprocesses.
        :param include_likelihood: Whether or not to add likelihood noise to
                                   the predicted underlying latent function f.
        :type include_likelihood: bool

        :returns: (mean, var):
            mean: posterior mean, a Numpy array, Nnew x self.input_dim
            var: posterior variance, a Numpy array, Nnew x 1 if full_cov=False,
                 Nnew x Nnew otherwise

            If full_cov and self.input_dim > 1, the return shape of var is
            Nnew x Nnew x self.input_dim. If self.input_dim == 1, the return
            shape is Nnew x Nnew. This is to allow for different normalizations
            of the output dimensions.

        Note: If you want the predictive quantiles (e.g. 95% confidence
        interval) use :py:func:"~GPy.core.gp.GP.predict_quantiles".