docs: def ec_matrix_vector(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 via successive summation of the following
      matrix vector products, p_1^t=p_0^t T,...,p_n^t=P_{n-1}^t T.
  
      Such a direct approach can become prohibetively expensive for
      large n.  In this case the sum can be computed more efficiently
      using an eigenvalue decomposition of T and applying the closed
      form expression for a finite geometric series to each of the
      eigenvalues.
  
      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:
          """Probability vector after (k=0) propagations"""
          p_k = 1.0 * p0
          """Sum of vectors after (k=0) propagations"""
          p_sum = 1.0 * p_k
          for k in range(n - 1):
              """Propagate one step p_{k} -> p_{k+1}"""
              p_k = np.dot(p_k, T)
              """Update sum"""
              p_sum += p_k
          """Expected counts"""
          EC = p_sum[:, np.newaxis] * T
          return EC