docs: def ec_geometric_series(p0, T, n):
      r"""Compute expected transition counts for Markov chain after n
      steps.
  
      Expected counts are computed according to ..math::
  
      E[C_{ij}^{(n)}]=\sum_{k=0}^{n-1} (p_0^t T^{k})_{i} p_{ij}
  
      The sum is computed using the eigenvalue decomposition of T and
      applying the expression for a finite geometric series to each of
      the eigenvalues.
  
      For small n the computation of the eigenvalue decomposition can be
      much more expensive than a direct computation. In this case it is
      beneficial to compute the expected counts using successively
      computed matrix vector products p_1^t=p_0^t T, ... as increments.
  
      Parameters
      ----------
      p0 : (M,) ndarray
          Starting (probability) vector of the chain.
      T : (M, M) ndarray
          Transition matrix of the chain.
      n : int
          Number of steps to take from initial state.
  
      Returns
      --------
      EC : (M, M) ndarray
          Expected value for transition counts after N steps.
  
      """
      if (n <= 0):
          EC = np.zeros(T.shape)
          return EC
      else:
          R, D, L = rdl_decomposition(T)
          w = np.diagonal(D)
          L = np.transpose(L)
  
          D_sum = np.diag(geometric_series(w, n - 1))
          T_sum = np.dot(np.dot(R, D_sum), np.conjugate(np.transpose(L)))
          p_sum = np.dot(p0, T_sum)
          EC = p_sum[:, np.newaxis] * T
          """Truncate imginary part - which is zero, but we want real
          return values"""
          EC = EC.real
          return EC