calculate the sensitivity matrix for index of the backward committor from A to B given transition matrix T.
    Parameters
    ----------
    T : numpy.ndarray shape = (n, n)
        Transition matrix
    A : array like
        List of integer state labels for set A
    B : array like
        List of integer state labels for set B
    index : entry of the committor 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.

def backward_committor_sensitivity(T, A, B, index):
    """
    calculate the sensitivity matrix for index of the backward committor from A to B given transition matrix T.
    Parameters
    ----------
    T : numpy.ndarray shape = (n, n)
        Transition matrix
    A : array like
        List of integer state labels for set A
    B : array like
        List of integer state labels for set B
    index : entry of the committor 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.
    """

    # This is really ugly to compute. The problem is, that changes in T induce changes in
    # the stationary distribution and so we need to add this influence, too
    # I implemented something which is correct, but don't ask me about the derivation

    n = len(T)

    trT = numpy.transpose(T)

    one = numpy.ones(n)
    eq = stationary_distribution(T)

    mEQ = numpy.diag(eq)
    mIEQ = numpy.diag(1.0 / eq)
    mSEQ = numpy.diag(1.0 / eq / eq)

    backT = numpy.dot(mIEQ, numpy.dot(trT, mEQ))

    qMat = forward_committor_sensitivity(backT, A, B, index)

    matA = trT - numpy.identity(n)
    matA = numpy.concatenate((matA, [one]))

    phiM = numpy.linalg.pinv(matA)

    phiM = phiM[:, 0:n]

    trQMat = numpy.transpose(qMat)

    d1 = numpy.dot(mSEQ, numpy.diagonal(numpy.dot(numpy.dot(trT, mEQ), trQMat), 0))
    d2 = numpy.diagonal(numpy.dot(numpy.dot(trQMat, mIEQ), trT), 0)

    psi1 = numpy.dot(d1, phiM)
    psi2 = numpy.dot(-d2, phiM)

    v1 = psi1 - one * numpy.dot(psi1, eq)
    v3 = psi2 - one * numpy.dot(psi2, eq)

    part1 = numpy.outer(eq, v1)
    part2 = numpy.dot(numpy.dot(mEQ, trQMat), mIEQ)
    part3 = numpy.outer(eq, v3)

    sensitivity = part1 + part2 + part3

    return sensitivity