code: def compute_upper_limit(mu_in, post, alpha=0.9):
      
      if 0 < alpha < 1:
          dp = integral_element(mu_in, post)
          high_idx = bisect.bisect_left(dp.cumsum() / dp.sum(), alpha)
          # if alpha is in (0,1] and post is non-negative, bisect_left
          # will always return an index in the range of mu since
          # post.cumsum()/post.sum() will always begin at 0 and end at 1
          mu_high = mu_in[high_idx]
      elif alpha == 1:
          mu_high = numpy.max(mu_in[post > 0])
      else:
          raise ValueError("Confidence level must be in (0,1].")
  
      return mu_high
docs: Returns the upper limit mu_high of confidence level alpha for a
      posterior distribution post on the given parameter mu.
      The posterior need not be normalized.