Group ``m`` (Markov) chains whose common :py:func:`.r_value` is
    less than ``critical_r`` in each of the D dimensions.

    :param means:

        (m x D) Matrix-like array; the mean value estimates.

    :param variances:

        (m x D) Matrix-like array; the variance estimates.

def r_group(means, variances, n, critical_r=2., approx=False):
    '''Group ``m`` (Markov) chains whose common :py:func:`.r_value` is
    less than ``critical_r`` in each of the D dimensions.

    :param means:

        (m x D) Matrix-like array; the mean value estimates.

    :param variances:

        (m x D) Matrix-like array; the variance estimates.

    '''
    assert len(means) == len(variances), \
    'Number of ``means`` (%i) does not match number of ``variances`` (%i)' % (len(means), len(variances))
    means = _np.asarray(means)
    variances  = _np.asarray(variances)
    assert means.ndim == 2, '``means`` must be matrix-like'
    assert variances.ndim == 2, '``variances`` must be 2-dimensional'
    assert means.shape[1] == variances.shape[1], \
    'Dimensionality of ``means`` (%i) and ``variances`` (%i) does not match' % (means.shape[1], variances.shape[1])

    groups = []

    for i in range(len(means)):
        assigned = False
        # try to assign component i to an existing group
        for group in groups:
            rows = group + [i]
            # R values for each parameter
            r_values = _np.array([r_value(means[rows, j], variances[rows, j], n, approx) for j in range(means.shape[1])])
            if _np.all(r_values < critical_r):
                # add to group if R value small enough
                group.append(i)
                assigned = True
                break
        # if component i has not been added to an existing group case create a new group
        if not assigned:
            groups.append([i])

    return groups