docs: def ep_pairwise(n_items, data, alpha, model="logit", max_iter=100,
          initial_state=None):
      """Compute a distribution of model parameters using the EP algorithm.
  
      This function computes an approximate Bayesian posterior probability
      distribution over model parameters, given pairwise-comparison data (see
      :ref:`data-pairwise`). It uses the expectation propagation algorithm, as
      presented, e.g., in [CG05]_.
  
      The prior distribution is assumed to be isotropic Gaussian with variance
      ``1 / alpha``. The posterior is approximated by a a general multivariate
      Gaussian distribution, described by a mean vector and a covariance matrix.
  
      Two different observation models are available. ``logit`` (default) assumes
      that pairwise-comparison outcomes follow from a Bradley-Terry model.
      ``probit`` assumes that the outcomes follow from Thurstone's model.
  
      Parameters
      ----------
      n_items : int
          Number of distinct items.
      data : list of lists
          Pairwise-comparison data.
      alpha : float
          Inverse variance of the (isotropic) prior.
      model : str, optional
          Observation model. Either "logit" or "probit".
      max_iter : int, optional
          Maximum number of iterations allowed.
      initial_state : tuple of array_like, optional
          Natural parameters used to initialize the EP algorithm.
  
      Returns
      -------
      mean : numpy.ndarray
          The mean vector of the approximate Gaussian posterior.
      cov : numpy.ndarray
          The covariance matrix of the approximate Gaussian posterior.
  
      Raises
      ------
      ValueError
          If the observation model is not "logit" or "probit".
      """
      if model == "logit":
          match_moments = _match_moments_logit
      elif model == "probit":
          match_moments = _match_moments_probit
      else:
          raise ValueError("unknown model '{}'".format(model))
      return _ep_pairwise(
              n_items, data, alpha, match_moments, max_iter, initial_state)