Pareto smoothed importance sampling (PSIS).

    Parameters
    ----------
    log_weights : array
        Array of size (n_samples, n_observations)
    reff : float
        relative MCMC efficiency, `ess / n`

    Returns
    -------
    lw_out : array
        Smoothed log weights
    kss : array
        Pareto tail indices

def psislw(log_weights, reff=1.0):
    """
    Pareto smoothed importance sampling (PSIS).

    Parameters
    ----------
    log_weights : array
        Array of size (n_samples, n_observations)
    reff : float
        relative MCMC efficiency, `ess / n`

    Returns
    -------
    lw_out : array
        Smoothed log weights
    kss : array
        Pareto tail indices
    """
    rows, cols = log_weights.shape

    log_weights_out = np.copy(log_weights, order="F")
    kss = np.empty(cols)

    # precalculate constants
    cutoff_ind = -int(np.ceil(min(rows / 5.0, 3 * (rows / reff) ** 0.5))) - 1
    cutoffmin = np.log(np.finfo(float).tiny)  # pylint: disable=no-member, assignment-from-no-return
    k_min = 1.0 / 3

    # loop over sets of log weights
    for i, x in enumerate(log_weights_out.T):
        # 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)  # pylint: disable=unbalanced-tuple-unpacking
        x_tail = x[tailinds]
        tail_len = len(x_tail)
        if tail_len <= 4:
            # not enough tail samples for gpdfit
            k = np.inf
        else:
            # order of tail samples
            x_tail_si = np.argsort(x_tail)
            # fit generalized Pareto distribution to the right tail samples
            x_tail = np.exp(x_tail) - expxcutoff
            k, sigma = _gpdfit(x_tail[x_tail_si])

            if k >= k_min:
                # no smoothing if short tail or GPD fit failed
                # compute ordered statistic for the fit
                sti = np.arange(0.5, tail_len) / tail_len
                smoothed_tail = _gpinv(sti, k, sigma)
                smoothed_tail = np.log(  # pylint: disable=assignment-from-no-return
                    smoothed_tail + expxcutoff
                )
                # place the smoothed tail into the output array
                x[tailinds[x_tail_si]] = smoothed_tail
                # truncate smoothed values to the largest raw weight 0
                x[x > 0] = 0
        # renormalize weights
        x -= _logsumexp(x)
        # store tail index k
        kss[i] = k

    return log_weights_out, kss