def predict(self, y, t=None, return_cov=True, return_var=False):
        """
        Compute the conditional predictive distribution of the model

        You must call :func:`GP.compute` before this method.

        Args:
            y (array[n]): The observations at coordinates ``x`` from
                :func:`GP.compute`.
            t (Optional[array[ntest]]): The independent coordinates where the
                prediction should be made. If this is omitted the coordinates
                will be assumed to be ``x`` from :func:`GP.compute` and an
                efficient method will be used to compute the prediction.
            return_cov (Optional[bool]): If ``True``, the full covariance
                matrix is computed and returned. Otherwise, only the mean
                prediction is computed. (default: ``True``)
            return_var (Optional[bool]): If ``True``, only return the diagonal
                of the predictive covariance; this will be faster to compute
                than the full covariance matrix. This overrides ``return_cov``
                so, if both are set to ``True``, only the diagonal is computed.
                (default: ``False``)

        Returns:
            ``mu``, ``(mu, cov)``, or ``(mu, var)`` depending on the values of
            ``return_cov`` and ``return_var``. These output values are:
            (a) **mu** ``(ntest,)``: mean of the predictive distribution,
            (b) **cov** ``(ntest, ntest)``: the predictive covariance matrix,
            and
            (c) **var** ``(ntest,)``: the diagonal elements of ``cov``.

        Raises:
            ValueError: For mismatched dimensions.

        """
        y = self._process_input(y)
        if len(y.shape) > 1:
            raise ValueError("dimension mismatch")

        if t is None:
            xs = self._t
        else:
            xs = np.ascontiguousarray(t, dtype=float)
            if len(xs.shape) > 1:
                raise ValueError("dimension mismatch")

        # Make sure that the model is computed
        self._recompute()

        # Compute the predictive mean.
        resid = y - self.mean.get_value(self._t)

        if t is None:
            alpha = self.solver.solve(resid).flatten()
            alpha = resid - (self._yerr**2 + self.kernel.jitter) * alpha
        elif not len(self._A):
            alpha = self.solver.predict(resid, xs)
        else:
            Kxs = self.get_matrix(xs, self._t)
            alpha = np.dot(Kxs, alpha)

        mu = self.mean.get_value(xs) + alpha
        if not (return_var or return_cov):
            return mu

        # Predictive variance.
        Kxs = self.get_matrix(xs, self._t)
        KxsT = np.ascontiguousarray(Kxs.T, dtype=np.float64)
        if return_var:
            var = -np.sum(KxsT*self.apply_inverse(KxsT), axis=0)
            var += self.kernel.get_value(0.0)
            return mu, var

        # Predictive covariance
        cov = self.kernel.get_value(xs[:, None] - xs[None, :])
        cov -= np.dot(Kxs, self.apply_inverse(KxsT))
        return mu, cov
Compute the conditional predictive distribution of the model

        You must call :func:`GP.compute` before this method.

        Args:
            y (array[n]): The observations at coordinates ``x`` from
                :func:`GP.compute`.
            t (Optional[array[ntest]]): The independent coordinates where the
                prediction should be made. If this is omitted the coordinates
                will be assumed to be ``x`` from :func:`GP.compute` and an
                efficient method will be used to compute the prediction.
            return_cov (Optional[bool]): If ``True``, the full covariance
                matrix is computed and returned. Otherwise, only the mean
                prediction is computed. (default: ``True``)
            return_var (Optional[bool]): If ``True``, only return the diagonal
                of the predictive covariance; this will be faster to compute
                than the full covariance matrix. This overrides ``return_cov``
                so, if both are set to ``True``, only the diagonal is computed.
                (default: ``False``)

        Returns:
            ``mu``, ``(mu, cov)``, or ``(mu, var)`` depending on the values of
            ``return_cov`` and ``return_var``. These output values are:
            (a) **mu** ``(ntest,)``: mean of the predictive distribution,
            (b) **cov** ``(ntest, ntest)``: the predictive covariance matrix,
            and
            (c) **var** ``(ntest,)``: the diagonal elements of ``cov``.

        Raises:
            ValueError: For mismatched dimensions.