r"""Calculate the sensitivity matrix for entry j the stationary
    distribution vector given transition matrix T.

    Parameters
    ----------
    T : numpy.ndarray shape = (n, n)
        Transition matrix
    j : int
        entry of stationary distribution for which the sensitivity is to be computed

    Returns
    -------
    x : ndarray, shape=(n, n)
        Sensitivity matrix for entry index around transition matrix T. Reversibility is not assumed.

    Remark
    ------
    Note, that this function uses a different normalization convention for the sensitivity compared to
    eigenvector_sensitivity. See there for further information.

def stationary_distribution_sensitivity(T, j):
    r"""Calculate the sensitivity matrix for entry j the stationary
    distribution vector given transition matrix T.

    Parameters
    ----------
    T : numpy.ndarray shape = (n, n)
        Transition matrix
    j : int
        entry of stationary distribution for which the sensitivity is to be computed

    Returns
    -------
    x : ndarray, shape=(n, n)
        Sensitivity matrix for entry index around transition matrix T. Reversibility is not assumed.

    Remark
    ------
    Note, that this function uses a different normalization convention for the sensitivity compared to
    eigenvector_sensitivity. See there for further information.
    """

    n = len(T)

    lEV = numpy.ones(n)
    rEV = stationary_distribution(T)
    eVal = 1.0

    T = numpy.transpose(T)

    vecA = numpy.zeros(n)
    vecA[j] = 1.0

    matA = T - eVal * numpy.identity(n)
    # normalize s.t. sum is one using rEV which is constant
    matA = numpy.concatenate((matA, [lEV]))

    phi = numpy.linalg.lstsq(numpy.transpose(matA), vecA, rcond=-1)
    phi = numpy.delete(phi[0], -1)

    sensitivity = -numpy.outer(rEV, phi) + numpy.dot(phi, rEV) * numpy.outer(rEV, lEV)

    return sensitivity