docs: def wavg(datalist, fast=False, prior=None, **fitterargs):
      """ Weighted average of |GVar|\s or arrays/dicts of |GVar|\s.
  
      The weighted average of ``N`` |GVar|\s ::
  
          xavg = wavg([g1, g2 ... gN])
  
      is what one obtains from a weighted least-squares fit of the
      collection of |GVar|\s to the one-parameter fit function ::
  
          def f(p):
              return N * [p]
  
      The average is the best-fit value for fit parameter ``p``.  |GVar|\s
      with smaller standard deviations carry more weight than those with
      larger standard deviations; and the averages take account of
      correlations between the |GVar|\s.
  
      ``wavg`` also works when each ``gi`` is an array of |GVar|\s or a
      dictionary whose values are |GVar|\s or arrays of |GVar|\s.
      Corresponding arrays in different ``gi``\s must have the same dimension,
      but can have different shapes (the overlapping components are
      averaged).  When the ``gi`` are dictionaries, they need not all have
      the same keys.
  
      Weighted averages can become costly when the number of random samples
      being averaged is large (100s or more). In such cases it might be useful
      to set parameter ``fast=True``. This causes ``wavg`` to estimate the
      weighted average by incorporating the random samples one at a time into a
      running average::
  
          result = datalist[0]
          for di in datalist[1:]:
              result = wavg([result, di], ...)
  
      This method can be much faster when ``len(datalist)`` is large, and gives
      the exact result when there are no correlations between different elements
      of list ``datalist``. The results are approximately correct when
      ``datalist[i]`` and ``datalist[j]`` are correlated for ``i!=j``.
  
      Args:
          datalist (list): The |GVar|\s to be averaged. ``datalist`` is
              a one-dimensional sequence of |GVar|\s, or of arrays of |GVar|\s,
              or of dictionaries containing |GVar|\s  and/or arrays of |GVar|\s.
              Corresponding arrays in different ``datalist[i]``\s must have the
              same dimension.
          fast (bool): If ``fast=True``, ``wavg`` averages the ``datalist[i]``
              sequentially. This can be much faster when averaging a large
              number of sampes but is only approximate if the different
              elements of ``datalist`` are correlated. Default is ``False``.
          fitterargs (dict): Additional arguments (e.g., ``svdcut``) for the
              :class:`lsqfit.nonlinear_fit` fitter used to do the averaging.
  
      Returns:
          The weighted average is returned as a |GVar| or an array of
          |GVar|\s or a dictionary of |GVar|\s and arrays of |GVar|\s.
          Results have the following extra attributes:
  
          **chi2** - ``chi**2`` for weighted average.
  
          **dof** - Effective number of degrees of freedom.
  
          **Q** - Quality factor `Q` (or *p-value*) for fit:
              the probability that the ``chi**2`` could have been larger,
              by chance, assuming that the data are all Gaussian and consistent
              with each other. Values smaller than 0.1 or so suggest that the
              data are not Gaussian or are inconsistent with each other. Also
              called the *p-value*.
  
          **time** - Time required to do average.
  
          **svdcorrection** - The *svd* corrections made to the data
              when ``svdcut`` is not ``None``.
  
          **fit** - Fit returned by :class:`lsqfit.nonlinear_fit`.
      """
      if prior is not None:
          datalist = list(datalist) + [prior]
          warnings.warn(
              'use of prior in lsqfit.wavg is deprecated',
              DeprecationWarning
              )
      if len(datalist) <= 0:
          return None
      elif len(datalist) == 1:
          if hasattr(datalist[0], 'keys'):
              return BufferDictWAvg(datalist[0], None)
          if numpy.shape(datalist[0]) == ():
              return GVarWAvg(datalist[0], None)
          else:
              return ArrayWAvg(numpy.asarray(datalist[0]), None)
      if fast:
          chi2 = dof = time = svdcorrection = 0
          ans = datalist[0]
          for i, di in enumerate(datalist[1:]):
              ans = wavg([ans, di], fast=False, **fitterargs)
              chi2 += ans.chi2
              dof += ans.dof
              time += ans.time
              svdcorrection += ans.svdcorrection
          ans.fit.dof = dof
          ans.fit.Q = _gammaQ(dof / 2., chi2 / 2.)
          ans.fit.chi2 = chi2
          ans.fit.time = time
          ans.fit.svdcorrection = svdcorrection
          return ans
      if hasattr(datalist[0], 'keys'):
          datashape = None
      else:
          datashape = numpy.shape(datalist[0])
          datalist = [{None:di} for di in datalist]
      # repack as a single dictionary
      p0shape = {}
      p0index = {}
      data = gvar.BufferDict()
      for i, di in enumerate(datalist):
          for k in di:
              data[k, i] = di[k]
              shape = numpy.shape(di[k])
              p0index[k, i] = tuple(slice(0, j) for j in shape)
              if k not in p0shape:
                  p0shape[k] = shape
              elif p0shape[k] != shape:
                  p0shape[k] = tuple(
                      max(j1, j2) for j1, j2 in zip(shape, p0shape[k])
                      )
      # calculate p0
      p0 = gvar.BufferDict()
      p0count = {}
      for k, i in data:
          if k not in p0:
              p0[k] = numpy.zeros(p0shape[k], float)
              p0count[k] = numpy.zeros(p0shape[k], float)
          if p0index[k, i] == ():
              p0[k] += data[k, i].mean
              p0count[k] += 1
          else:
              p0[k][p0index[k, i]] += gvar.mean(data[k, i])
              p0count[k][p0index[k, i]] += 1.
      for k in p0:
          p0[k] /= p0count[k]
      # set up fit
      def fcn(p):
          ans = gvar.BufferDict()
          for k, i in data:
              shape = data[k, i].shape
              if shape == ():
                  ans[k, i] = p[k]
              else:
                  ans[k, i] = p[k][p0index[k, i]]
          return ans
      fit = lsqfit.nonlinear_fit(data=data, fcn=fcn, p0=p0, **fitterargs)
      if datashape is None:
          return BufferDictWAvg(fit.p, fit)
      elif datashape == ():
          return GVarWAvg(fit.p[None], fit)
      else:
          return ArrayWAvg(fit.p[None], fit)