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.

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