docs: def _calc_CI(
      sampler,
      modelidx=0,
      confs=[3, 1],
      last_step=False,
      e_range=None,
      e_npoints=100,
      threads=None,
  ):
      """Calculate confidence interval.
      """
      from scipy import stats
  
      # If we are computing the samples for the confidence intervals, we need at
      # least one sample to constrain the highest confidence band
      # 1 sigma -> 6 samples
      # 2 sigma -> 43 samples
      # 3 sigma -> 740 samples
      # 4 sigma -> 31574 samples
      # 5 sigma -> 3488555 samples
      # We limit it to 1000 samples and warn that it might not be enough
      if e_range:
          maxconf = np.max(confs)
          minsamples = min(100, int(1 / stats.norm.cdf(-maxconf) + 1))
          if minsamples > 1000:
              log.warning(
                  "In order to sample the confidence band for {0} sigma,"
                  " {1} new samples need to be computed, but we are limiting"
                  " it to 1000 samples, so the confidence band might not be"
                  " well constrained."
                  " Consider reducing the maximum"
                  " confidence significance or using the samples stored in"
                  " the sampler by setting e_range"
                  " to None".format(maxconf, minsamples)
              )
              minsamples = 1000
      else:
          minsamples = None
  
      modelx, model = _read_or_calc_samples(
          sampler,
          modelidx,
          last_step=last_step,
          e_range=e_range,
          e_npoints=e_npoints,
          n_samples=minsamples,
          threads=threads,
      )
  
      nwalkers = len(model) - 1
      CI = []
      for conf in confs:
          fmin = stats.norm.cdf(-conf)
          fmax = stats.norm.cdf(conf)
          ymin, ymax = [], []
          for fr, y in ((fmin, ymin), (fmax, ymax)):
              nf = int((fr * nwalkers))
              for i in six.moves.range(len(modelx)):
                  ysort = np.sort(model[:, i])
                  y.append(ysort[nf])
  
          # create an array from lists ymin and ymax preserving units
          CI.append((u.Quantity(ymin), u.Quantity(ymax)))
  
      return modelx, CI