docs: def _ep_pairwise(
          n_items, comparisons, alpha, match_moments, max_iter, initial_state):
      """Compute a distribution of model parameters using the EP algorithm.
  
      Raises
      ------
      RuntimeError
          If the algorithm does not converge after ``max_iter`` iterations.
      """
      # Static variable that allows to check the # of iterations after the call.
      _ep_pairwise.iterations = 0
      m = len(comparisons)
      prior_inv = alpha * np.eye(n_items)
      if initial_state is None:
          # Initially, mean and covariance come from the prior.
          mean = np.zeros(n_items)
          cov = (1 / alpha) * np.eye(n_items)
          # Initialize the natural params in the function space.
          tau = np.zeros(m)
          nu = np.zeros(m)
          # Initialize the natural params in the space of thetas.
          prec = np.zeros((n_items, n_items))
          xs = np.zeros(n_items)
      else:
          tau, nu = initial_state
          mean, cov, xs, prec = _init_ws(
                  n_items, comparisons, prior_inv, tau, nu)
      for _ in range(max_iter):
          _ep_pairwise.iterations += 1
          # Keep a copy of the old parameters for convergence testing.
          tau_old = np.array(tau, copy=True)
          nu_old = np.array(nu, copy=True)
          for i in nprand.permutation(m):
              a, b = comparisons[i]
              # Update mean and variance in function space.
              f_var = cov[a,a] + cov[b,b] - 2 * cov[a,b]
              f_mean = mean[a] - mean[b]
              # Cavity distribution.
              tau_tot = 1.0 / f_var
              nu_tot = tau_tot * f_mean
              tau_cav = tau_tot - tau[i]
              nu_cav = nu_tot - nu[i]
              cov_cav = 1.0 / tau_cav
              mean_cav = cov_cav * nu_cav
              # Moment matching.
              logpart, dlogpart, d2logpart = match_moments(mean_cav, cov_cav)
              # Update factor params in the function space.
              tau[i] = -d2logpart / (1 + d2logpart / tau_cav)
              delta_tau = tau[i] - tau_old[i]
              nu[i] = ((dlogpart - (nu_cav / tau_cav) * d2logpart)
                       / (1 + d2logpart / tau_cav))
              delta_nu = nu[i] - nu_old[i]
              # Update factor params in the weight space.
              prec[(a, a, b, b), (a, b, a, b)] += delta_tau * MAT_ONE_FLAT
              xs[a] += delta_nu
              xs[b] -= delta_nu
              # Update mean and covariance.
              if abs(delta_tau) > 0:
                  phi = -1.0 / ((1.0 / delta_tau) + f_var) * MAT_ONE
                  upd_mat = cov.take([a, b], axis=0)
                  cov = cov + upd_mat.T.dot(phi).dot(upd_mat)
              mean = cov.dot(xs)
          # Recompute the global parameters for stability.
          cov = inv_posdef(prior_inv + prec)
          mean = cov.dot(xs)
          if _converged((tau, nu), (tau_old, nu_old)):
              return mean, cov
      raise RuntimeError(
              "EP did not converge after {} iterations".format(max_iter))