code: def moments_match_ep(self, Y_i, tau_i, v_i, Y_metadata_i=None):
          
          if Y_i == 1:
              sign = 1.
          elif Y_i == 0 or Y_i == -1:
              sign = -1
          else:
              raise ValueError("bad value for Bernoulli observation (0, 1)")
          if isinstance(self.gp_link, link_functions.Probit):
              z = sign*v_i/np.sqrt(tau_i**2 + tau_i)
              phi_div_Phi = derivLogCdfNormal(z)
              log_Z_hat = logCdfNormal(z)
  
              mu_hat = v_i/tau_i + sign*phi_div_Phi/np.sqrt(tau_i**2 + tau_i)
              sigma2_hat = 1./tau_i - (phi_div_Phi/(tau_i**2+tau_i))*(z+phi_div_Phi)
  
          elif isinstance(self.gp_link, link_functions.Heaviside):
              z = sign*v_i/np.sqrt(tau_i)
              phi_div_Phi = derivLogCdfNormal(z)
              log_Z_hat = logCdfNormal(z)
              mu_hat = v_i/tau_i + sign*phi_div_Phi/np.sqrt(tau_i)
              sigma2_hat = (1. - a*phi_div_Phi - np.square(phi_div_Phi))/tau_i
          else:
              #TODO: do we want to revert to numerical quadrature here?
              raise ValueError("Exact moment matching not available for link {}".format(self.gp_link.__name__))
  
          # TODO: Output log_Z_hat instead of Z_hat (needs to be change in all others likelihoods)
          return np.exp(log_Z_hat), mu_hat, sigma2_hat
docs: Moments match of the marginal approximation in EP algorithm
  
          :param i: number of observation (int)
          :param tau_i: precision of the cavity distribution (float)
          :param v_i: mean/variance of the cavity distribution (float)