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.

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