ed by corners `a` and `b`.

If a subclass does not override this method, then a default error estimator is
used. This estimates the error as ``|est - refined_est|`` where ``est`` is
``estimate(f, a, b)`` and ``refined_est`` is the sum of
``estimate(f, a_k, b_k)`` where ``a_k, b_k`` are the coordinates of each
subregion of the region described by ``a`` and ``b``. In the 1D case, this
is equivalent to comparing the integral over an entire interval ``[a, b]`` to
the sum of the integrals over the left and right subintervals, ``[a, (a+b)/2]``
and ``[(a+b)/2, b]``.

Parameters
----------
f : callable
    Function to estimate error for. `f` must have the signature::
        f(x : ndarray, \*args) -> ndarray

    `f` should accept arrays `x` of shape::
        (npoints, ndim)

    and output arrays of shape::
        (npoints, output_dim_1, ..., output_dim_n)

    In this case, `estimate` will return arrays of shape::
        (output_dim_1, ..., output_dim_n)
a, b : ndarray
    Lower and upper limits of integration as rank-1 arrays specifying the left
    and right endpoints of the intervals being integrated over. Infinite limits
    are currently not supported.
args : tuple, optional
    Additional positional args passed to `f`, if any.

Returns
-------
err_est : ndarray
    Result of error estimation. If `f` returns arrays of shape
    ``(npoints, output_dim_1, ..., output_dim_n)``, then `est` will be
    of shape ``(output_dim_1, ..., output_dim_n)``.
r