code: def bayesianUpperLimit(self, alpha, steps=1.e5, plot=False):
          
          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)
          pdf = np.exp(f(x) / 2.)
          cut = (pdf / np.max(pdf)) > 1.e-10
          x = x[cut]
          pdf = pdf[cut]
          #pdf /= pdf[0]
          #forbidden = np.nonzero(pdf < 1.e-10)[0]
          #if len(forbidden) > 0:
          #    index = forbidden[0] # Numeric stability
          #    x = x[0: index]
          #    pdf = pdf[0: index]
          cdf = np.cumsum(pdf)
          cdf /= cdf[-1]
          cdf_reflect = scipy.interpolate.interp1d(cdf, x)
  
          return cdf_reflect(alpha)
docs: Compute one-sided upper limit using Bayesian Method of Helene.
          Several methods of increasing numerical stability have been implemented.