docs: def psislw(lw, Reff=1.0, overwrite_lw=False):
      """Pareto smoothed importance sampling (PSIS).
  
      Parameters
      ----------
      lw : ndarray
          Array of size n x m containing m sets of n log weights. It is also
          possible to provide one dimensional array of length n.
  
      Reff : scalar, optional
          relative MCMC efficiency ``N_eff / N``
  
      overwrite_lw : bool, optional
          If True, the input array `lw` is smoothed in-place, assuming the array
          is F-contiguous. By default, a new array is allocated.
  
      Returns
      -------
      lw_out : ndarray
          smoothed log weights
      kss : ndarray
          Pareto tail indices
  
      """
      if lw.ndim == 2:
          n, m = lw.shape
      elif lw.ndim == 1:
          n = len(lw)
          m = 1
      else:
          raise ValueError("Argument `lw` must be 1 or 2 dimensional.")
      if n <= 1:
          raise ValueError("More than one log-weight needed.")
  
      if overwrite_lw and lw.flags.f_contiguous:
          # in-place operation
          lw_out = lw
      else:
          # allocate new array for output
          lw_out = np.copy(lw, order='F')
  
      # allocate output array for kss
      kss = np.empty(m)
  
      # precalculate constants
      cutoff_ind = - int(np.ceil(min(0.2 * n, 3 * np.sqrt(n / Reff)))) - 1
      cutoffmin = np.log(np.finfo(float).tiny)
      logn = np.log(n)
      k_min = 1/3
  
      # loop over sets of log weights
      for i, x in enumerate(lw_out.T if lw_out.ndim == 2 else lw_out[None, :]):
          # improve numerical accuracy
          x -= np.max(x)
          # sort the array
          x_sort_ind = np.argsort(x)
          # divide log weights into body and right tail
          xcutoff = max(
              x[x_sort_ind[cutoff_ind]],
              cutoffmin
          )
          expxcutoff = np.exp(xcutoff)
          tailinds, = np.where(x > xcutoff)
          x2 = x[tailinds]
          n2 = len(x2)
          if n2 <= 4:
              # not enough tail samples for gpdfitnew
              k = np.inf
          else:
              # order of tail samples
              x2si = np.argsort(x2)
              # fit generalized Pareto distribution to the right tail samples
              np.exp(x2, out=x2)
              x2 -= expxcutoff
              k, sigma = gpdfitnew(x2, sort=x2si)
          if k >= k_min and not np.isinf(k):
              # no smoothing if short tail or GPD fit failed
              # compute ordered statistic for the fit
              sti = np.arange(0.5, n2)
              sti /= n2
              qq = gpinv(sti, k, sigma)
              qq += expxcutoff
              np.log(qq, out=qq)
              # place the smoothed tail into the output array
              x[tailinds[x2si]] = qq
              # truncate smoothed values to the largest raw weight 0
              x[x > 0] = 0
          # renormalize weights
          x -= sumlogs(x)
          # store tail index k
          kss[i] = k
  
      # If the provided input array is one dimensional, return kss as scalar.
      if lw_out.ndim == 1:
          kss = kss[0]
  
      return lw_out, kss