docs: def bootstrapped_fit_iter(self, n=None, datalist=None, **kargs):
          """ Iterator that returns bootstrap copies of a fit.
  
          A bootstrap analysis involves three steps: 1) make a large number
          of "bootstrap copies" of the original input data and prior that differ
          from each other by random amounts characteristic of the underlying
          randomness in the original data; 2) repeat the entire fit analysis
          for each bootstrap copy of the data, extracting fit results from
          each; and 3) use the variation of the fit results from bootstrap
          copy to bootstrap copy to determine an approximate probability
          distribution (possibly non-gaussian) for the fit parameters and/or
          functions of them: the results from each bootstrap fit are samples
          from that distribution.
  
          Bootstrap copies of the data for step 2 are provided in
          ``datalist``. If ``datalist`` is ``None``, they are generated
          instead from the means and covariance matrix of the fit data
          (assuming gaussian statistics). The maximum number of bootstrap
          copies considered is specified by ``n`` (``None`` implies no
          limit).
  
          Variations in the best-fit parameters (or functions of them)
          from bootstrap fit to bootstrap fit define the probability
          distributions for those quantities. For example, one could use the
          following code to analyze the distribution of function ``g(p)``
          of the fit parameters::
  
              fit = nonlinear_fit(...)
  
              ...
  
              glist = []
              for bsfit in fit.bootstrapped_fit_iter(
                  n=100, datalist=datalist,
                  ):
                  glist.append(g(bsfit.pmean))
  
              ... analyze samples glist[i] from g(p) distribution ...
  
          This code generates ``n=100`` samples ``glist[i]`` from the
          probability distribution of ``g(p)``. If everything is Gaussian,
          the mean and standard deviation of ``glist[i]`` should agree
          with ``g(fit.p).mean`` and ``g(fit.p).sdev``.
  
          Args:
              n (int): Maximum number of iterations if ``n`` is not ``None``;
                  otherwise there is no maximum.
              datalist (iter): Collection of bootstrap ``data`` sets for fitter.
              kargs (dict): Overrides arguments in original fit.
  
          Returns:
              Iterator that returns an |nonlinear_fit| object
              containing results from the fit to the next data set in
              ``datalist``.
          """
          fargs = dict(fitter=self.fitter, fcn=self.fcn)
          fargs.update(self.fitterargs)
          fargs['p0'] = self.pmean
          fargs['prior'] = self.prior
          for k in kargs:
              fargs[k] = kargs[k]
          prior = fargs['prior']
          del fargs['prior']
          if datalist is None:
              x = self.x
              y = self.y
              if prior is None:
                  for yb in _gvar.bootstrap_iter(y, n):
                      fit = nonlinear_fit(data=(x, yb), prior=None, **fargs)
                      yield fit
              else:
                  g = _gvar.BufferDict(y=y.flat, prior=prior.flat)
                  for gb in _gvar.bootstrap_iter(g, n):
                      yb = _reformat(y, buf=gb['y'])
                      priorb = _reformat(prior, buf=gb['prior'])
                      fit = nonlinear_fit(data=(x, yb), prior=priorb, **fargs)
                      yield fit
          else:
              if prior is None:
                  i = 0
                  for datab in datalist:
                      i += 1
                      if n is not None and i > n:
                          break
                      fit = nonlinear_fit(data=datab, prior=None, **fargs)
                      yield fit
              else:
                  piter = _gvar.bootstrap_iter(prior)
                  i = 0
                  for datab in datalist:
                      i += 1
                      if n is not None and i > n:
                          break
                      fit = nonlinear_fit(data=datab, prior=next(piter), **fargs)
                      yield fit