docs: def r_value(means, variances, n, approx=False):
      '''Calculate the Gelman-Rubin R value (Chapter 2.2 in [GR92]_).
  
      The R value can be used to quantify mixing of "multiple iterative
      simulations" (e.g. Markov Chains) in parameter space.  An R value
      "close to one" indicates that all chains explored the same region
      of the parameter.
  
      .. note::
  
          The R value is defined only in *one* dimension.
  
      :param means:
  
          Vector-like array; the sample mean of each chain.
  
      :param variances:
  
          Vector-like array; the sample variance of each chain.
  
      '''
      # use same variable names as in [GR92]
      # means is \bar{x}_i
      # variances is s_i^2
  
      means     = _np.asarray(means)
      variances = _np.asarray(variances)
  
      assert means.ndim == 1, '``means`` must be vector-like'
      assert variances.ndim == 1, '``variances`` must be vector-like'
      assert len(means) == len(variances), \
      'Number of ``means`` (%i) does not match number of ``variances`` (%i)' %( len(means), len(variances) )
  
      m = len(means)
  
      x_bar    = _np.average(means)
      B_over_n = ((means - x_bar)**2).sum() / (m - 1)
  
      W = _np.average(variances)
  
      # var_estimate is \hat{\sigma}^2
      var_estimate = (n - 1) / n  *  W + B_over_n
  
      if approx:
          return var_estimate / W
  
      V = var_estimate + B_over_n / m
  
      # calculate the three terms of \hat{var}(\hat{V}) (equation (4) in [GR92]
      # var_V is \hat{var}(\hat{V})
      tmp_cov_matrix = _np.cov(variances, means)
      # third term
      var_V = _np.cov(variances, means**2)[1,0] - 2. * x_bar * tmp_cov_matrix[1,0]
      var_V *= 2. * (m + 1) * (n - 1) / (m * m * n)
      # second term (merged n in denominator into ``B_over_n``)
      var_V += ((m + 1) / m)**2 * 2. / (m - 1) * B_over_n * B_over_n
      # first term
      var_V += ((n - 1) / n)**2 / m * tmp_cov_matrix[0,0]
  
      df = 2. * V**2 / var_V
  
      if df <= 2.:
          return _np.inf
  
      return V / W * df / (df - 2)