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
def psislw ( log_weights , reff = 1.0 ) : rows , cols = log_weights . shape log_weights_out = np . copy ( log_weights , order = "F" ) kss = np . empty ( cols ) cutoff_ind = - int ( np . ceil ( min ( rows / 5.0 , 3 * ( rows / reff ) ** 0.5 ) ) ) - 1 cutoffmin = np . log ( np . finfo ( float ) . tiny ) k_min = 1.0 / 3 for i , x in enumerate ( log_weights_out . T ) : x -= np . max ( x ) x_sort_ind = np . argsort ( x ) xcutoff = max ( x [ x_sort_ind [ cutoff_ind ] ] , cutoffmin ) expxcutoff = np . exp ( xcutoff ) tailinds , = np . where ( x > xcutoff ) x_tail = x [ tailinds ] tail_len = len ( x_tail ) if tail_len <= 4 : k = np . inf else : x_tail_si = np . argsort ( x_tail ) x_tail = np . exp ( x_tail ) - expxcutoff k , sigma = _gpdfit ( x_tail [ x_tail_si ] ) if k >= k_min : sti = np . arange ( 0.5 , tail_len ) / tail_len smoothed_tail = _gpinv ( sti , k , sigma ) smoothed_tail = np . log ( smoothed_tail + expxcutoff ) x [ tailinds [ x_tail_si ] ] = smoothed_tail x [ x > 0 ] = 0 x -= _logsumexp ( x ) kss [ i ] = k return log_weights_out , kss
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