docs: def bootstrap_histogram_2D(
          xvalues, yvalues, xbins, ybins,
          boundaries=[None, None], xsigma=None, ysigma=None,
          normalisation=False, number_bootstraps=None):
      '''
      Calculates a 2D histogram of data, allowing for normalisation and
      bootstrap sampling
  
      :param numpy.ndarray xvalues:
          Data values of the first variable
  
      :param numpy.ndarray yvalues:
          Data values of the second variable
  
      :param numpy.ndarray xbins:
          Bin edges for the first variable
  
      :param numpy.ndarray ybins:
          Bin edges for the second variable
  
      :param list boundaries:
          List of (Lower, Upper) tuples corresponding to the bounds of the
          two data sets
  
      :param numpy.ndarray xsigma:
          Error values (standard deviatons) on first variable
  
      :param numpy.ndarray ysigma:
          Error values (standard deviatons) on second variable
  
      :param bool normalisation:
          If True then returns the histogram as a density function
  
      :param int number_bootstraps:
          Number of bootstraps
  
      :param returns:
          2-D histogram of data
      '''
      if (xsigma is None and ysigma is None) or not number_bootstraps:
          # No sampling - return simple 2-D histrogram
          #output = np.histogram2d(xvalues, yvalues, bins=[xbins, ybins])[0]
          output = hmtk_histogram_2D(xvalues, yvalues, bins=(xbins, ybins))
          if normalisation:
              output = output / float(np.sum(output))
          return output
  
      else:
          if xsigma is None:
              xsigma = np.zeros(len(xvalues), dtype=float)
          if ysigma is None:
              ysigma = np.zeros(len(yvalues), dtype=float)
          temp_hist = np.zeros(
              [len(xbins) - 1, len(ybins) - 1, number_bootstraps],
              dtype=float)
          for iloc in range(0, number_bootstraps):
              xsample = sample_truncated_gaussian_vector(xvalues, xsigma,
                                                         boundaries[0])
              ysample = sample_truncated_gaussian_vector(yvalues, ysigma,
                                                         boundaries[0])
  
              # temp_hist[:, :, iloc] = np.histogram2d(xsample,
              #                                       ysample,
              #                                       bins=[xbins, ybins])[0]
              temp_hist[:, :, iloc] = hmtk_histogram_2D(xsample,
                                                        ysample,
                                                        bins=(xbins, ybins))
          if normalisation:
              output = np.sum(temp_hist, axis=2)
              output = output / np.sum(output)
          else:
              output = np.sum(temp_hist, axis=2) / float(number_bootstraps)
          return output