docs: 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