docs: def confidenceInterval(self, alpha=0.6827, steps=1.e5, plot=False):
          """
          Compute two-sided confidence interval by taking x-values corresponding to the largest PDF-values first.
          """
          x_dense, y_dense = self.densify()
          y_dense -= np.max(y_dense) # Numeric stability
          f = scipy.interpolate.interp1d(x_dense, y_dense, kind='linear')
          x = np.linspace(0., np.max(x_dense), steps)
          # ADW: Why does this start at 0, which often outside the input range?
          # Wouldn't starting at xmin be better:
          #x = np.linspace(np.min(x_dense), np.max(x_dense), steps)
          pdf = np.exp(f(x) / 2.)
          cut = (pdf / np.max(pdf)) > 1.e-10
          x = x[cut]
          pdf = pdf[cut]
  
          sorted_pdf_indices = np.argsort(pdf)[::-1] # Indices of PDF in descending value
          cdf = np.cumsum(pdf[sorted_pdf_indices])
          cdf /= cdf[-1]
          sorted_pdf_index_max = np.argmin((cdf - alpha)**2)
          x_select = x[sorted_pdf_indices[0: sorted_pdf_index_max]]
              
          return np.min(x_select), np.max(x_select)