def evaluate_likelihood(self, y, f, hyp=None):
        """
        Evaluate the Gaussian CDF likelihood model,
            Φ(yₙfₙ) = (1 + erf(yₙfₙ / √2)) / 2
        for erf(z) = (2/√π) ∫ exp(-x²) dx, where the integral is over [0, z]
        Can be used to evaluate Q quadrature points when performing moment matching.
        :param y: observed data yₙ ϵ {-1, +1} [scalar]
        :param f: latent function value fₙ [Q, 1]
        :param hyp: dummy input, Probit has no hyperparameters
        :return:
            Φ(yₙfₙ) [Q, 1]
        """
        return (1.0 + erf(y * f / np.sqrt(2.0))) / 2.0  # Φ(z)