def conditional_moments(self, f, hyp=None):
        """
        The first two conditional moments of a Probit likelihood are:
            E[yₙ|fₙ] = Φ(fₙ)
            Var[yₙ|fₙ] = Φ(fₙ) (1 - Φ(fₙ))
            where Φ(fₙ) = (1 + erf(fₙ / √2)) / 2
        """
        # TODO: not working
        # phi = (1.0 + erf(f / np.sqrt(2.0))) / 2.0
        # phi = self.link_fn(f)
        # phi = erfc(f / np.sqrt(2.0)) - 1.0
        phi = self.evaluate_likelihood(1.0, f)
        return phi, phi * (1.0 - phi)