docs: def compute_threshold(smoothed_errors, error_buffer, sd_limit=12.0):
      """Helper method for `extract_anomalies` method.
      Calculates the epsilon (threshold) for anomalies.
      """
      mu = np.mean(smoothed_errors)
      sigma = np.std(smoothed_errors)
  
      max_epsilon = 0
      sd_threshold = sd_limit
  
      # The treshold is determined dynamically by testing multiple Zs.
      # z is drawn from an ordered set of positive values representing the
      # number of standard deviations above mean(smoothed_errors)
  
      # here we iterate in increments of 0.5 on the range that the NASA paper found to be good
      for z in np.arange(2.5, sd_limit, 0.5):
          epsilon = mu + (sigma * z)
          below_epsilon, below_indices, above_epsilon = [], [], []
  
          for i in range(len(smoothed_errors)):
              e = smoothed_errors[i]
              if e < epsilon:
                  # save to compute delta mean and delta std
                  # these are important for epsilon calculation
                  below_epsilon.append(e)
                  below_indices.append(i)
  
              if e > epsilon:
                  # above_epsilon values are anomalies
                  for j in range(0, error_buffer):
                      if (i + j) not in above_epsilon and (i + j) < len(smoothed_errors):
                          above_epsilon.append(i + j)
  
                      if (i - j) not in above_epsilon and (i - j) >= 0:
                          above_epsilon.append(i - j)
  
          if len(above_epsilon) == 0:
              continue
  
          # generate sequences
          above_epsilon = sorted(list(set(above_epsilon)))
          groups = [list(group) for group in mit.consecutive_groups(above_epsilon)]
          above_sequences = [(g[0], g[-1]) for g in groups if not g[0] == g[-1]]
  
          mean_perc_decrease = (mu - np.mean(below_epsilon)) / mu
          sd_perc_decrease = (sigma - np.std(below_epsilon)) / sigma
          epsilon = (mean_perc_decrease + sd_perc_decrease) /\
                    (len(above_sequences)**2 + len(above_epsilon))
  
          # update the largest epsilon we've seen so far
          if epsilon > max_epsilon:
              sd_threshold = z
              max_epsilon = epsilon
  
      # sd_threshold can be multiplied by sigma to get epsilon
      return max_epsilon, sd_threshold