Classify the 2D ``data`` according to the ``responsibility`` and
    make a scatter plot of each data point with the color of the
    component it is most likely from. The ``responsibility`` is
    normalized internally such that each row sums to unity.

    :param data:

        matrix-like; one row = one 2D sample

    :param responsibility:

        matrix-like; one row = probabilities that sample n is from
        1st, 2nd, ... component. The number of rows has to agree with ``data``

    :param cmap:

        colormap; defines how component indices are mapped to the
        color of the data points

def plot_responsibility(data, responsibility,
                        cmap='nipy_spectral'):
    '''Classify the 2D ``data`` according to the ``responsibility`` and
    make a scatter plot of each data point with the color of the
    component it is most likely from. The ``responsibility`` is
    normalized internally such that each row sums to unity.

    :param data:

        matrix-like; one row = one 2D sample

    :param responsibility:

        matrix-like; one row = probabilities that sample n is from
        1st, 2nd, ... component. The number of rows has to agree with ``data``

    :param cmap:

        colormap; defines how component indices are mapped to the
        color of the data points

    '''
    import numpy as np
    from matplotlib import pyplot as plt
    from matplotlib.cm import get_cmap

    data = np.asarray(data)
    responsibility = np.asarray(responsibility)

    assert data.ndim == 2
    assert responsibility.ndim == 2

    D = data.shape[1]
    N = data.shape[0]
    K = responsibility.shape[1]

    assert D == 2
    assert N == responsibility.shape[0]

    # normalize responsibility so each row sums to one
    inv_row_sum = 1.0 / np.einsum('nk->n', responsibility)
    responsibility = np.einsum('n,nk->nk', inv_row_sum, responsibility)

    # index of the most likely component for each sample
    indicators = np.argmax(responsibility, axis=1)

    # same color range as in plot_mixture
    if K > 1:
        point_colors = indicators / (K - 1) * _max_color
    else:
        point_colors = np.zeros(N)
    plt.scatter(data.T[0], data.T[1], c=point_colors, cmap=cmap)