docs: def chi2(g1, g2=None, svdcut=1e-12, nocorr=False):
      """ Compute chi**2 of ``g1-g2``.
  
      ``chi**2`` is a measure of whether the multi-dimensional
      Gaussian distributions ``g1`` and ``g2`` (dictionaries or arrays)
      agree with each other --- that is, do their means agree
      within errors for corresponding elements. The probability is high
      if ``chi2(g1,g2)/chi2.dof`` is of order 1 or smaller.
  
      Usually ``g1`` and ``g2`` are dictionaries with the same keys,
      where ``g1[k]`` and ``g2[k]`` are |GVar|\s or arrays of
      |GVar|\s having the same shape. Alternatively ``g1`` and ``g2``
      can be |GVar|\s, or arrays of |GVar|\s having the same shape.
  
      One of ``g1`` or ``g2`` can contain numbers instead of |GVar|\s,
      in which case ``chi**2`` is a measure of the likelihood that
      the numbers came from the distribution specified by the other
      argument.
  
      One or the other of ``g1`` or ``g2`` can be missing keys, or missing
      elements from arrays. Only the parts of ``g1`` and ``g2`` that
      overlap are used. Also setting ``g2=None`` is equivalent to replacing its
      elements by zeros.
  
      ``chi**2`` is computed from the inverse of the covariance matrix
      of ``g1-g2``. The matrix inversion can be sensitive to roundoff
      errors. In such cases, SVD cuts can be applied by setting
      parameters ``svdcut``; see the documentation
      for :func:`gvar.svd`, which is used to apply the cut.
  
      The return value is the ``chi**2``. Extra attributes attached to this
      value give additional information:
  
      - **dof** --- Number of degrees of freedom (that is, the number of variables
        compared).
  
      - **Q** --- The probability that the ``chi**2`` could have been larger,
        by chance, even if ``g1`` and ``g2`` agree. Values smaller than 0.1
        or so suggest that they do not agree. Also called the *p-value*.
      """
      # customized class for answer
      class ans(float):
          def __new__(cls, chi2, dof, Q):
              return float.__new__(cls, chi2)
          def __init__(self, chi2, dof, Q):
              self.dof = dof
              self.Q = Q
              self.chi2 = chi2
  
      # leaving nocorr (turn off correlations) undocumented because I
      #   suspect I will remove it
      if g2 is None:
          diff = BufferDict(g1).buf if hasattr(g1, 'keys') else numpy.asarray(g1).flatten()
      elif hasattr(g1, 'keys') and hasattr(g2, 'keys'):
          # g1 and g2 are dictionaries
          g1 = BufferDict(g1)
          g2 = BufferDict(g2)
          diff = BufferDict()
          keys = set(g1.keys())
          keys = keys.intersection(g2.keys())
          for k in keys:
              g1k = g1[k]
              g2k = g2[k]
              shape = tuple(
                  [min(s1,s2) for s1, s2 in zip(numpy.shape(g1k), numpy.shape(g2k))]
                  )
              diff[k] = numpy.zeros(shape, object)
              if len(shape) == 0:
                  diff[k] = g1k - g2k
              else:
                  for i in numpy.ndindex(shape):
                      diff[k][i] = g1k[i] - g2k[i]
          diff = diff.buf
      elif not hasattr(g1, 'keys') and not hasattr(g2, 'keys'):
          # g1 and g2 are arrays or scalars
          g1 = numpy.asarray(g1)
          g2 = numpy.asarray(g2)
          shape = tuple(
              [min(s1,s2) for s1, s2 in zip(numpy.shape(g1), numpy.shape(g2))]
              )
          diff = numpy.zeros(shape, object)
          if len(shape) == 0:
              diff = numpy.array(g1 - g2)
          else:
              for i in numpy.ndindex(shape):
                  diff[i] = g1[i] - g2[i]
          diff = diff.flatten()
      else:
          # g1 and g2 are something else
          raise ValueError(
              'cannot compute chi**2 for types ' + str(type(g1)) + ' ' +
              str(type(g2))
              )
      dof = diff.size
      if dof == 0:
          return ans(0.0, 0, 0)
      if nocorr:
          # ignore correlations
          chi2 = numpy.sum(mean(diff) ** 2 / var(diff))
          dof = len(diff)
      else:
          diffmod, i_wgts = svd(diff, svdcut=svdcut, wgts=-1)
          diffmean = mean(diffmod)
          i, wgts = i_wgts[0]
          chi2 = 0.0
          if len(i) > 0:
              chi2 += numpy.sum((diffmean[i] * wgts) ** 2)
          for i, wgts in i_wgts[1:]:
              chi2 += numpy.sum(wgts.dot(diffmean[i]) ** 2)
          dof = sum(len(wgts) for i, wgts in i_wgts)
      Q = gammaQ(dof/2., chi2/2.)
      return ans(chi2, dof=dof, Q=Q)