docs: def gpdfitnew(x, sort=True, sort_in_place=False, return_quadrature=False):
      """Estimate the paramaters for the Generalized Pareto Distribution (GPD)
  
      Returns empirical Bayes estimate for the parameters of the two-parameter
      generalized Parato distribution given the data.
  
      Parameters
      ----------
      x : ndarray
          One dimensional data array
  
      sort : bool or ndarray, optional
          If known in advance, one can provide an array of indices that would
          sort the input array `x`. If the input array is already sorted, provide
          False. If True (default behaviour), the array is sorted internally.
  
      sort_in_place : bool, optional
          If `sort` is True and `sort_in_place` is True, the array is sorted
          in-place (False by default).
  
      return_quadrature : bool, optional
          If True, quadrature points and weight `ks` and `w` of the marginal posterior distribution of k are also calculated and returned. False by
          default.
  
      Returns
      -------
      k, sigma : float
          estimated parameter values
  
      ks, w : ndarray
          Quadrature points and weights of the marginal posterior distribution
          of `k`. Returned only if `return_quadrature` is True.
  
      Notes
      -----
      This function returns a negative of Zhang and Stephens's k, because it is
      more common parameterisation.
  
      """
      if x.ndim != 1 or len(x) <= 1:
          raise ValueError("Invalid input array.")
  
      # check if x should be sorted
      if sort is True:
          if sort_in_place:
              x.sort()
              xsorted = True
          else:
              sort = np.argsort(x)
              xsorted = False
      elif sort is False:
          xsorted = True
      else:
          xsorted = False
  
      n = len(x)
      PRIOR = 3
      m = 30 + int(np.sqrt(n))
  
      bs = np.arange(1, m + 1, dtype=float)
      bs -= 0.5
      np.divide(m, bs, out=bs)
      np.sqrt(bs, out=bs)
      np.subtract(1, bs, out=bs)
      if xsorted:
          bs /= PRIOR * x[int(n/4 + 0.5) - 1]
          bs += 1 / x[-1]
      else:
          bs /= PRIOR * x[sort[int(n/4 + 0.5) - 1]]
          bs += 1 / x[sort[-1]]
  
      ks = np.negative(bs)
      temp = ks[:,None] * x
      np.log1p(temp, out=temp)
      np.mean(temp, axis=1, out=ks)
  
      L = bs / ks
      np.negative(L, out=L)
      np.log(L, out=L)
      L -= ks
      L -= 1
      L *= n
  
      temp = L - L[:,None]
      np.exp(temp, out=temp)
      w = np.sum(temp, axis=1)
      np.divide(1, w, out=w)
  
      # remove negligible weights
      dii = w >= 10 * np.finfo(float).eps
      if not np.all(dii):
          w = w[dii]
          bs = bs[dii]
      # normalise w
      w /= w.sum()
  
      # posterior mean for b
      b = np.sum(bs * w)
      # Estimate for k, note that we return a negative of Zhang and
      # Stephens's k, because it is more common parameterisation.
      temp = (-b) * x # pylint: disable=invalid-unary-operand-type
      np.log1p(temp, out=temp)
      k = np.mean(temp)
      if return_quadrature:
          np.negative(x, out=temp)
          temp = bs[:, None] * temp
          np.log1p(temp, out=temp)
          ks = np.mean(temp, axis=1)
      # estimate for sigma
      sigma = -k / b * n / (n - 0)
      # weakly informative prior for k
      a = 10
      k = k * n / (n+a) + a * 0.5 / (n+a)
      if return_quadrature:
          ks *= n / (n+a)
          ks += a * 0.5 / (n+a)
  
      if return_quadrature:
          return k, sigma, ks, w
      else:
          return k, sigma