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)