code: def _gpdfit(x):
      
      prior_bs = 3
      prior_k = 10
      len_x = len(x)
      m_est = 30 + int(len_x ** 0.5)
  
      b_ary = 1 - np.sqrt(m_est / (np.arange(1, m_est + 1, dtype=float) - 0.5))
      b_ary /= prior_bs * x[int(len_x / 4 + 0.5) - 1]
      b_ary += 1 / x[-1]
  
      k_ary = np.log1p(-b_ary[:, None] * x).mean(axis=1)  # pylint: disable=no-member
      len_scale = len_x * (np.log(-(b_ary / k_ary)) - k_ary - 1)
      weights = 1 / np.exp(len_scale - len_scale[:, None]).sum(axis=1)
  
      # remove negligible weights
      real_idxs = weights >= 10 * np.finfo(float).eps
      if not np.all(real_idxs):
          weights = weights[real_idxs]
          b_ary = b_ary[real_idxs]
      # normalise weights
      weights /= weights.sum()
  
      # posterior mean for b
      b_post = np.sum(b_ary * weights)
      # estimate for k
      k_post = np.log1p(-b_post * x).mean()  # pylint: disable=invalid-unary-operand-type,no-member
      # add prior for k_post
      k_post = (len_x * k_post + prior_k * 0.5) / (len_x + prior_k)
      sigma = -k_post / b_post
  
      return k_post, sigma
docs: Estimate the parameters for the Generalized Pareto Distribution (GPD).
  
      Empirical Bayes estimate for the parameters of the generalized Pareto
      distribution given the data.
  
      Parameters
      ----------
      x : array
          sorted 1D data array
  
      Returns
      -------
      k : float
          estimated shape parameter
      sigma : float
          estimated scale parameter