Add convenience function for computing expectation-specific diagnostics

[sethaxen]

The PSIS paper notes it is sometimes useful to compute h-specific (function/expectation-specific) diagnostics, which more specifically diagnose the quality of importance sampling for the expectation in question. From discussions with @avehtari, this can be done with the following procedure. Given importance ratios r(θ) and function evaluations h(θ), we:

1. Compute w(θ) as r(θ) with a Pareto-smoothed upper tail.
2. Compute v(θ) = h(θ) r(θ). Fit the GPD to upper tails of both v(θ) and -v(θ) to get kₕᵤ and kₕₗ. Compute kₕ=max(kₕᵤ, kₕₗ).
3. Warn for either k > 0.7 or kₕ > 0.7
4. Return w(θ), k, and kₕ. The user would then estimate 𝔼[h] as sum(h .* w) ./ sum(w).

Likewise, we can return h-specific ESS estimates by simply replacing w in the ESS computation with v.

A few things to consider:

1. Providing this requires decoupling diagnosing from smoothing. i.e. for some uses one might just want to fit the GPD and check the diagnostic without modifying any weights. (Edit: see #26 26)
2. PSISResult currently can't accommodate 3 smoothing results. Perhaps it should be extended, or a different object should be returned if h is provided.
3. Currently the API assumes users sample from the proposal distribution once and then smooth for many different target distributions simultaneously. While this is what is needed for LOO, I suspect that diagnosing for many different h's simultaneously could be more useful. Perhaps the API should support both.
4. RE (2), since Monte Carlo is just a special case of IS where all weights are uniform, the same function could be used to diagnose all univariate marginals of a posterior sample to identify problems with estimating means and variances of parameters. This should be supported.
5. the moment-matching paper (https://arxiv.org/abs/1906.08850, see #25) uses a similar formulation with v. Between this feature and that, some changes may need to be made to the internals.