diff --git a/Project.toml b/Project.toml
index 078f7720..83e65190 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,24 +1,23 @@
 name = "PSIS"
 uuid = "ce719bf2-d5d0-4fb9-925d-10a81b42ad04"
-version = "0.9.8"
+version = "0.10.0-DEV"
 authors = ["Seth Axen <seth.axen@gmail.com> and contributors"]
 
 [deps]
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 DimensionalData = "0.24.2, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30"
 Distributions = "0.25.81"
+DocStringExtensions = "0.9.2"
 JLD2 = "0.4.46, 0.5, 0.6"
 LinearAlgebra = "1"
 LogExpFunctions = "0.3.3"
-Plots = "1.10.1"
 Printf = "1"
-RecipesBase = "1"
 ReferenceTests = "0.9, 0.10"
 StableRNGs = "1"
 Statistics = "1"
@@ -29,10 +28,9 @@ DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DimensionalData", "Distributions", "Logging", "JLD2", "Plots", "ReferenceTests", "StableRNGs", "Test"]
+test = ["DimensionalData", "Distributions", "Logging", "JLD2", "ReferenceTests", "StableRNGs", "Test"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 8b322147..27224cdf 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -5,7 +5,6 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 MCMCDiagnosticTools = "be115224-59cd-429b-ad48-344e309966f0"
 PSIS = "ce719bf2-d5d0-4fb9-925d-10a81b42ad04"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
@@ -14,4 +13,3 @@ Documenter = "1"
 DocumenterCitations = "1.2.1"
 DocumenterInterLinks = "1"
 MCMCDiagnosticTools = "0.3.2"
-Plots = "1.10.1"
diff --git a/docs/make.jl b/docs/make.jl
index 3366675e..07aa3653 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,4 +1,4 @@
-using Distributions, Plots, PSIS
+using Distributions, PSIS
 using Documenter
 using DocumenterCitations
 using DocumenterInterLinks
@@ -22,7 +22,6 @@ makedocs(;
     ),
     pages=[
         "Home" => "index.md",
-        "Plotting" => "plotting.md",
         "API" => "api.md",
         "Internal" => "internal.md",
         "References" => "references.md",
diff --git a/docs/src/api.md b/docs/src/api.md
index 5e7b69c5..107de24f 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -5,12 +5,4 @@
 ```@docs
 PSISResult
 psis
-ess_is
-```
-
-## Plotting
-
-```@docs
-PSISPlots
-PSISPlots.paretoshapeplot
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
index 5cc5b4f7..953e0eb7 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -32,6 +32,3 @@ result # hide
 ```
 
 As indicated by the warnings, this is a poor choice of a proposal distribution, and estimates are unlikely to converge (see [`PSISResult`](@ref) for an explanation of the shape thresholds).
-
-When running PSIS with many parameters, it is useful to plot the Pareto shape values to diagnose convergence.
-See [Plotting PSIS results](@ref) for examples.
diff --git a/docs/src/plotting.md b/docs/src/plotting.md
deleted file mode 100644
index 9c391d83..00000000
--- a/docs/src/plotting.md
+++ /dev/null
@@ -1,28 +0,0 @@
-# Plotting PSIS results
-
-PSIS.jl includes plotting recipes for [`PSISResult`](@ref) using any Plots.jl backend, as well as the utility plotting function [`PSISPlots.paretoshapeplot`](@ref).
-
-We demonstrate this with a simple example.
-
-```@example 1
-using Random # hide
-Random.seed!(42) # hide
-using PSIS, Distributions
-proposal = Normal()
-target = TDist(7)
-ndraws, nchains, nparams = (1_000, 1, 20)
-x = rand(proposal, ndraws, nchains, nparams)
-log_ratios = logpdf.(target, x) .- logpdf.(proposal, x)
-result = psis(log_ratios)
-```
-
-## Plots.jl
-
-`PSISResult` objects can be plotted directly:
-
-```@example 1
-using Plots
-plot(result; showlines=true, marker=:+, legend=false, linewidth=2)
-```
-
-This is equivalent to calling `PSISPlots.paretoshapeplot(result; kwargs...)`.
diff --git a/src/PSIS.jl b/src/PSIS.jl
index f249a936..67d8601e 100644
--- a/src/PSIS.jl
+++ b/src/PSIS.jl
@@ -1,17 +1,16 @@
 module PSIS
 
+using DocStringExtensions: FIELDS
 using LogExpFunctions: LogExpFunctions
 using Printf: @sprintf
 using Statistics: Statistics
 
-export PSISPlots
 export PSISResult
-export psis, psis!, ess_is
+export psis
 
 include("utils.jl")
 include("generalized_pareto.jl")
 include("core.jl")
 include("ess.jl")
-include("recipes/plots.jl")
 
 end
diff --git a/src/core.jl b/src/core.jl
index 7f2cc3d9..3adf2922 100644
--- a/src/core.jl
+++ b/src/core.jl
@@ -1,42 +1,26 @@
 # range, description, condition
-const SHAPE_DIAGNOSTIC_CATEGORIES = (
+const KHAT_DIAGNOSTIC_CATEGORIES = (
     ("(-Inf, 0.5]", "good", ≤(0.5)),
     ("(0.5, 0.7]", "okay", x -> 0.5 < x ≤ 0.7),
     ("(0.7, 1]", "bad", x -> 0.7 < x ≤ 1),
     ("(1, Inf)", "very bad", >(1)),
     ("——", "failed", isnan),
 )
-const BAD_SHAPE_SUMMARY = "Resulting importance sampling estimates are likely to be unstable."
-const VERY_BAD_SHAPE_SUMMARY = "Corresponding importance sampling estimates are likely to be unstable and are unlikely to converge with additional samples."
+const BAD_KHAT_SUMMARY = "Resulting importance sampling estimates are likely to be unstable."
+const VERY_BAD_KHAT_SUMMARY = "Corresponding importance sampling estimates are likely to be unstable and are unlikely to converge with additional samples."
 
 """
     PSISResult
 
 Result of Pareto-smoothed importance sampling (PSIS) using [`psis`](@ref).
 
-# Properties
-
-  - `log_weights`: un-normalized Pareto-smoothed log weights
-  - `weights`: normalized Pareto-smoothed weights (allocates a copy)
-  - `pareto_shape`: Pareto ``k=ξ`` shape parameter
-  - `nparams`: number of parameters in `log_weights`
-  - `ndraws`: number of draws in `log_weights`
-  - `nchains`: number of chains in `log_weights`
-  - `reff`: the ratio of the effective sample size of the unsmoothed importance ratios and
-    the actual sample size.
-  - `ess`: estimated effective sample size of estimate of mean using smoothed importance
-    samples (see [`ess_is`](@ref))
-  - `tail_length`: length of the upper tail of `log_weights` that was smoothed
-  - `tail_dist`: the generalized Pareto distribution that was fit to the tail of
-    `log_weights`. Note that the tail weights are scaled to have a maximum of 1, so
-    `tail_dist * exp(maximum(log_ratios))` is the corresponding fit directly to the tail of
-    `log_ratios`.
-  - `normalized::Bool`:indicates whether `log_weights` are log-normalized along the sample
-    dimensions.
+# Fields
+
+$(FIELDS)
 
 # Diagnostic
 
-The `pareto_shape` parameter ``k=ξ`` of the generalized Pareto distribution `tail_dist` can
+The `pareto_khat` parameter ``k=ξ`` of the generalized Pareto distribution can
 be used to diagnose reliability and convergence of estimates using the importance weights
 [VehtariSimpson2021](@citep).
 
@@ -53,74 +37,37 @@ be used to diagnose reliability and convergence of estimates using the importanc
     exists. The convergence rate is close to zero, and bias can be large with practical
     sample sizes.
 
-See [`PSISPlots.paretoshapeplot`](@ref) for a diagnostic plot.
-
 # References
 
   - [VehtariSimpson2021](@cite) Vehtari et al. JMLR 25:72 (2021).
 """
-struct PSISResult{T,W<:AbstractArray{T},R,L,D}
+struct PSISResult{T,W<:AbstractArray{T},K,E}
+    """Un-normalized Pareto-smoothed log importance weights."""
     log_weights::W
-    reff::R
-    tail_length::L
-    tail_dist::D
-    normalized::Bool
-end
-
-function Base.propertynames(r::PSISResult)
-    return [fieldnames(typeof(r))..., :weights, :nparams, :ndraws, :nchains, :pareto_shape]
-end
-
-function Base.getproperty(r::PSISResult, k::Symbol)
-    if k === :weights
-        log_weights = getfield(r, :log_weights)
-        getfield(r, :normalized) && return exp.(log_weights)
-        return LogExpFunctions.softmax(log_weights; dims=_sample_dims(log_weights))
-    elseif k === :nparams
-        log_weights = getfield(r, :log_weights)
-        return if ndims(log_weights) == 1
-            1
-        else
-            param_dims = _param_dims(log_weights)
-            prod(Base.Fix1(size, log_weights), param_dims; init=1)
-        end
-    elseif k === :ndraws
-        log_weights = getfield(r, :log_weights)
-        return size(log_weights, 1)
-    elseif k === :nchains
-        log_weights = getfield(r, :log_weights)
-        return size(log_weights, 2)
-    end
-    k === :pareto_shape && return pareto_shape(r)
-    k === :ess && return ess_is(r)
-    return getfield(r, k)
+    """Pareto k-hat (shape) diagnostic."""
+    pareto_khat::K
+    """Estimated effective sample size for the importance weights."""
+    ess_is::E
 end
 
 function Base.show(io::IO, ::MIME"text/plain", r::PSISResult)
-    npoints = r.nparams
-    nchains = r.nchains
-    println(
-        io, "PSISResult with $(r.ndraws) draws, $nchains chains, and $npoints parameters"
-    )
-    return _print_pareto_shape_summary(io, r; newline_at_end=false)
+    (ndraws, nchains, nparams) = _sample_param_sizes(r.log_weights)
+    println(io, "PSISResult with $ndraws draws, $nchains chains, and $nparams parameters")
+    return _print_pareto_khat_summary(io, r; newline_at_end=false)
 end
 
-function pareto_shape_summary(r::PSISResult; kwargs...)
-    return _print_pareto_shape_summary(stdout, r; kwargs...)
-end
-
-function _print_pareto_shape_summary(io::IO, r::PSISResult; kwargs...)
-    k = as_array(pareto_shape(r))
-    ess = as_array(ess_is(r))
-    npoints = r.nparams
-    rows = map(SHAPE_DIAGNOSTIC_CATEGORIES) do (range, desc, cond)
+function _print_pareto_khat_summary(io::IO, r::PSISResult; kwargs...)
+    k = as_array(pareto_khat(r))
+    ess_is = as_array(r.ess_is)
+    (_, _, nparams) = _sample_param_sizes(r.log_weights)
+    rows = map(KHAT_DIAGNOSTIC_CATEGORIES) do (range, desc, cond)
         inds = findall(cond, k)
         count = length(inds)
-        perc = 100 * count / npoints
+        perc = 100 * count / nparams
         ess_min = if count == 0 || desc == "failed"
-            oftype(first(ess), NaN)
+            oftype(first(ess_is), NaN)
         else
-            minimum(view(ess, inds))
+            minimum(view(ess_is, inds))
         end
         return (range=range, desc=desc, count_perc=(count, perc), ess_min=ess_min)
     end
@@ -143,7 +90,7 @@ function _print_pareto_shape_summary(io::IO, r::PSISResult; kwargs...)
         floor(Int, log10(maximum(r -> r.count_perc[2], rows))) + 6,
     ]
 
-    println(io, "Pareto shape (k) diagnostic values:")
+    println(io, "Pareto k-hat diagnostic summary:")
     printstyled(
         io,
         col_padding,
@@ -165,7 +112,12 @@ function _print_pareto_shape_summary(io::IO, r::PSISResult; kwargs...)
         format = formats[r.desc]
         printstyled(io, _pad_left(count, col_widths[3]); format...)
         printstyled(io, " ", _pad_right(perc_str, col_widths[4]); format...)
-        print(io, col_delim_tot, isfinite(r.ess_min) ? floor(Int, r.ess_min) : "——")
+        ess_min_string = if isfinite(r.ess_min) && r.desc ∈ ("good", "okay")
+            string(floor(Int, r.ess_min))
+        else
+            "——"
+        end
+        print(io, col_delim_tot, ess_min_string)
     end
     return nothing
 end
@@ -174,35 +126,30 @@ _pad_left(s, nchars) = " "^(nchars - length("$s")) * "$s"
 _pad_right(s, nchars) = "$s" * " "^(nchars - length("$s"))
 
 """
-    psis(log_ratios, reff = 1.0; kwargs...) -> PSISResult
-    psis!(log_ratios, reff = 1.0; kwargs...) -> PSISResult
+    psis(log_ratios, r_eff = 1.0; kwargs...) -> PSISResult
 
-Compute Pareto smoothed importance sampling (PSIS) log weights [VehtariSimpson2021](@citep).
-
-While `psis` computes smoothed log weights out-of-place, `psis!` smooths them in-place.
+Compute Pareto smoothed importance sampling (PSIS) log-weights [VehtariSimpson2021](@citep).
 
 # Arguments
 
   - `log_ratios`: an array of logarithms of importance ratios, with size
     `(draws, [chains, [parameters...]])`, where `chains>1` would be used when chains are
     generated using Markov chain Monte Carlo.
-  - `reff::Union{Real,AbstractArray}`: the ratio(s) of effective sample size of
-    `log_ratios` and the actual sample size `reff = ess/(draws * chains)`, used to account
+  - `r_eff::Union{Real,AbstractArray}`: the ratio(s) of effective sample size of
+    `log_ratios` and the actual sample size `r_eff = ess/(draws * chains)`, used to account
     for autocorrelation, e.g. due to Markov chain Monte Carlo. If an array, it must have the
     size `(parameters...,)` to match `log_ratios`.
 
 # Keywords
 
   - `warn=true`: If `true`, warning messages are delivered
-  - `normalize=true`: If `true`, the log-weights will be log-normalized so that
-    `exp.(log_weights)` sums to 1 along the sample dimensions.
 
 # Returns
 
   - `result`: a [`PSISResult`](@ref) object containing the results of the Pareto-smoothing.
 
-A warning is raised if the Pareto shape parameter ``k ≥ 0.7``. See [`PSISResult`](@ref) for
-details and [`PSISPlots.paretoshapeplot`](@ref) for a diagnostic plot.
+A warning is raised if the Pareto k-hat parameter ``k ≥ 0.7``. See [`PSISResult`](@ref) for
+details.
 
 # Examples
 
@@ -219,12 +166,12 @@ julia> x = rand(proposal, 1_000, 1, 30);  # (ndraws, nchains, nparams)
 julia> log_ratios = @. logpdf(target, x) - logpdf(proposal, x);
 
 julia> result = psis(log_ratios)
-┌ Warning: 9 parameters had Pareto shape values 0.7 < k ≤ 1. Resulting importance sampling estimates are likely to be unstable.
+┌ Warning: 9 parameters had Pareto k-hat values in (0.7, 1]. Resulting importance sampling estimates are likely to be unstable.
 └ @ PSIS ~/.julia/packages/PSIS/...
-┌ Warning: 1 parameters had Pareto shape values k > 1. Corresponding importance sampling estimates are likely to be unstable and are unlikely to converge with additional samples.
+┌ Warning: 1 parameters had Pareto k-hat > 1. Corresponding importance sampling estimates are likely to be unstable and are unlikely to converge with additional samples.
 └ @ PSIS ~/.julia/packages/PSIS/...
 PSISResult with 1000 draws, 1 chains, and 30 parameters
-Pareto shape (k) diagnostic values:
+Pareto k-hat diagnostic summary:
                         Count       Min. ESS
  (-Inf, 0.5]  good       7 (23.3%)  959
   (0.5, 0.7]  okay      13 (43.3%)  938
@@ -238,15 +185,15 @@ If the draws were generated using MCMC, we can compute the relative efficiency u
 ```jldoctest psis
 julia> using MCMCDiagnosticTools
 
-julia> reff = ess(log_ratios; kind=:basic, split_chains=1, relative=true);
+julia> r_eff = ess(log_ratios; kind=:basic, split_chains=1, relative=true);
 
-julia> result = psis(log_ratios, reff)
-┌ Warning: 9 parameters had Pareto shape values 0.7 < k ≤ 1. Resulting importance sampling estimates are likely to be unstable.
+julia> result = psis(log_ratios, r_eff)
+┌ Warning: 9 parameters had Pareto k-hat values in (0.7, 1]. Resulting importance sampling estimates are likely to be unstable.
 └ @ PSIS ~/.julia/packages/PSIS/...
-┌ Warning: 1 parameters had Pareto shape values k > 1. Corresponding importance sampling estimates are likely to be unstable and are unlikely to converge with additional samples.
+┌ Warning: 1 parameters had Pareto k-hat > 1. Corresponding importance sampling estimates are likely to be unstable and are unlikely to converge with additional samples.
 └ @ PSIS ~/.julia/packages/PSIS/...
 PSISResult with 1000 draws, 1 chains, and 30 parameters
-Pareto shape (k) diagnostic values:
+Pareto k-hat diagnostic summary:
                         Count       Min. ESS
  (-Inf, 0.5]  good       9 (30.0%)  806
   (0.5, 0.7]  okay      11 (36.7%)  842
@@ -258,30 +205,34 @@ Pareto shape (k) diagnostic values:
 
   - [VehtariSimpson2021](@cite) Vehtari et al. JMLR 25:72 (2021).
 """
-psis, psis!
+psis
 
-function psis(logr, reff=1; kwargs...)
+function psis(logr, r_eff=1; warn::Bool=true, kwargs...)
     T = float(eltype(logr))
     logw = similar(logr, T)
     copyto!(logw, logr)
-    return psis!(logw, reff; kwargs...)
+    _, khat = _psis!(logw, r_eff; warn, kwargs...)
+    ess = ess_is(logw, r_eff)
+    result = PSISResult(logw, khat, ess)
+    warn && check_pareto_khat(result)
+    return result
 end
 
-function psis!(logw::AbstractVecOrMat, reff=1; normalize::Bool=true, warn::Bool=true)
+function _psis!(logw::AbstractVecOrMat, r_eff; warn::Bool=true)
     T = typeof(float(one(eltype(logw))))
-    if length(reff) != 1
-        throw(DimensionMismatch("`reff` has length $(length(reff)) but must have length 1"))
+    if length(r_eff) != 1
+        throw(
+            DimensionMismatch("`r_eff` has length $(length(r_eff)) but must have length 1")
+        )
     end
-    warn && check_reff(reff)
+    warn && check_r_eff(r_eff)
     S = length(logw)
-    reff_val = first(reff)
-    M = tail_length(reff_val, S)
+    r_eff_val = first(r_eff)
+    M = tail_length(r_eff_val, S)
     if M < 5
         warn &&
             @warn "$M tail draws is insufficient to fit the generalized Pareto distribution. Total number of draws should in general exceed 25."
-        _maybe_log_normalize!(logw, normalize)
-        tail_dist_failed = GeneralizedPareto(0, T(NaN), T(NaN))
-        return PSISResult(logw, reff_val, M, tail_dist_failed, normalize)
+        return logw, T(NaN)
     end
     perm = partialsortperm(logw, (S - M):S)
     cutoff_ind = perm[1]
@@ -291,89 +242,75 @@ function psis!(logw::AbstractVecOrMat, reff=1; normalize::Bool=true, warn::Bool=
     if !all(isfinite, logw_tail)
         warn &&
             @warn "Tail contains non-finite values. Generalized Pareto distribution cannot be reliably fit."
-        _maybe_log_normalize!(logw, normalize)
-        tail_dist_failed = GeneralizedPareto(0, T(NaN), T(NaN))
-        return PSISResult(logw, reff_val, M, tail_dist_failed, normalize)
+        return logw, T(NaN)
     end
-    _, tail_dist = psis_tail!(logw_tail, logu)
-    warn && check_pareto_shape(tail_dist)
-    _maybe_log_normalize!(logw, normalize)
-    return PSISResult(logw, reff_val, M, tail_dist, normalize)
+    _, tail_dist = _psis_tail!(logw_tail, logu)
+    k = T(pareto_khat(tail_dist))
+    return logw, k
 end
-function psis!(logw::AbstractMatrix, reff=1; kwargs...)
-    result = psis!(vec(logw), reff; kwargs...)
+function _psis!(logw::AbstractMatrix, r_eff; kwargs...)
+    _, k = _psis!(vec(logw), r_eff; kwargs...)
     # unflatten log_weights
-    return PSISResult(
-        logw, result.reff, result.tail_length, result.tail_dist, result.normalized
-    )
+    return logw, k
 end
-function psis!(logw::AbstractArray, reff=1; normalize::Bool=true, warn::Bool=true)
+function _psis!(logw::AbstractArray, r_eff; warn::Bool=true)
     T = typeof(float(one(eltype(logw))))
     # if an array defines custom indices (e.g. AbstractDimArray), we preserve them
     param_axes = _param_axes(logw)
     param_shape = map(length, param_axes)
-    if !(length(reff) == 1 || size(reff) == param_shape)
+    if !(length(r_eff) == 1 || size(r_eff) == param_shape)
         throw(
             DimensionMismatch(
-                "`reff` has shape $(size(reff)) but must have same shape as the parameter axes $(param_shape)",
+                "`r_eff` has shape $(size(r_eff)) but must have same shape as the parameter axes $(param_shape)",
             ),
         )
     end
-    check_reff(reff)
+    check_r_eff(r_eff)
 
     # allocate containers
-    reffs = similar(logw, eltype(reff), param_axes)
-    reffs .= reff
-    tail_lengths = similar(logw, Int, param_axes)
-    tail_dists = similar(logw, GeneralizedPareto{T}, param_axes)
+    r_effs = similar(logw, eltype(r_eff), param_axes)
+    r_effs .= r_eff
+    khats = similar(logw, T, param_axes)
 
     # call psis! in parallel for all parameters
     Threads.@threads for i in _eachparamindex(logw)
         logw_i = _selectparam(logw, i)
-        result_i = psis!(logw_i, reffs[i]; normalize=normalize, warn=false)
-        tail_lengths[i] = result_i.tail_length
-        tail_dists[i] = result_i.tail_dist
+        _, k = _psis!(logw_i, r_effs[i]; warn=false)
+        khats[i] = k
     end
 
-    # combine results
-    result = PSISResult(logw, reffs, tail_lengths, map(identity, tail_dists), normalize)
-
-    # warn for bad shape
-    warn && check_pareto_shape(result)
-    return result
+    return logw, khats
 end
 
-pareto_shape(dist::GeneralizedPareto) = dist.k
-pareto_shape(r::PSISResult) = pareto_shape(getfield(r, :tail_dist))
-pareto_shape(dists) = map(pareto_shape, dists)
+pareto_khat(dist::GeneralizedPareto) = dist.k
+pareto_khat(r::PSISResult) = r.pareto_khat
 
-function check_reff(reff)
-    isvalid = all(reff) do r
+function check_r_eff(r_eff)
+    isvalid = all(r_eff) do r
         return isfinite(r) && r > 0
     end
-    isvalid || @warn "All values of `reff` should be finite, but some are not."
+    isvalid || @warn "All values of `r_eff` should be finite, but some are not."
     return nothing
 end
 
-check_pareto_shape(result::PSISResult) = check_pareto_shape(result.tail_dist)
-function check_pareto_shape(dist::GeneralizedPareto)
-    k = pareto_shape(dist)
+check_pareto_khat(result::PSISResult) = check_pareto_khat(result.pareto_khat)
+function check_pareto_khat(k::Real)
     if k > 1
-        @warn "Pareto shape k = $(@sprintf("%.2g", k)) > 1. $VERY_BAD_SHAPE_SUMMARY"
+        @warn "Pareto k-hat = $(@sprintf("%.2g", k)) > 1. $VERY_BAD_KHAT_SUMMARY"
     elseif k > 0.7
-        @warn "Pareto shape k = $(@sprintf("%.2g", k)) > 0.7. $BAD_SHAPE_SUMMARY"
+        @warn "Pareto k-hat = $(@sprintf("%.2g", k)) > 0.7. $BAD_KHAT_SUMMARY"
     end
     return nothing
 end
-function check_pareto_shape(dists::AbstractArray{<:GeneralizedPareto})
-    nnan = count(isnan ∘ pareto_shape, dists)
-    ngt07 = count(>(0.7) ∘ pareto_shape, dists)
-    ngt1 = iszero(ngt07) ? ngt07 : count(>(1) ∘ pareto_shape, dists)
+function check_pareto_khat(ks::AbstractArray{<:Real})
+    nnan = count(isnan, ks)
+    ngt07 = count(>(0.7), ks)
+    ngt1 = iszero(ngt07) ? ngt07 : count(>(1), ks)
     if ngt07 > ngt1
-        @warn "$(ngt07 - ngt1) parameters had Pareto shape values 0.7 < k ≤ 1. $BAD_SHAPE_SUMMARY"
+        @warn "$(ngt07 - ngt1) parameters had Pareto k-hat values in (0.7, 1]. $BAD_KHAT_SUMMARY"
     end
     if ngt1 > 0
-        @warn "$ngt1 parameters had Pareto shape values k > 1. $VERY_BAD_SHAPE_SUMMARY"
+        @warn "$ngt1 parameters had Pareto k-hat > 1. $VERY_BAD_KHAT_SUMMARY"
     end
     if nnan > 0
         @warn "For $nnan parameters, the generalized Pareto distribution could not be fit to the tail draws. Total number of draws should in general exceed 25, and the tail draws must be finite."
@@ -381,14 +318,14 @@ function check_pareto_shape(dists::AbstractArray{<:GeneralizedPareto})
     return nothing
 end
 
-function tail_length(reff, S)
+function tail_length(r_eff, S)
     max_length = cld(S, 5)
-    (isfinite(reff) && reff > 0) || return max_length
-    min_length = ceil(Int, 3 * sqrt(S / reff))
+    (isfinite(r_eff) && r_eff > 0) || return max_length
+    min_length = ceil(Int, 3 * sqrt(S / r_eff))
     return min(max_length, min_length)
 end
 
-function psis_tail!(logw, logμ)
+function _psis_tail!(logw, logμ)
     T = eltype(logw)
     logw_max = logw[end]
     # to improve numerical stability, we first shift the log-weights to have a maximum of 0,
@@ -397,7 +334,7 @@ function psis_tail!(logw, logμ)
     w_scaled = (logw .= exp.(logw .- logw_max) .- μ_scaled)
     tail_dist = fit_gpd(w_scaled; prior_adjusted=true, sorted=true)
     # undo the scaling
-    k = pareto_shape(tail_dist)
+    k = pareto_khat(tail_dist)
     if isfinite(k)
         p = uniform_probabilities(T, length(logw))
         @inbounds for i in eachindex(logw, p)
diff --git a/src/ess.jl b/src/ess.jl
index 9cc45416..40e60ee1 100644
--- a/src/ess.jl
+++ b/src/ess.jl
@@ -1,44 +1,27 @@
 """
-    ess_is(weights; reff=1)
+    ess_is(w; r_eff=1)
 
 Estimate effective sample size (ESS) for importance sampling over the sample dimensions.
 
-Given normalized weights ``w_{1:n}``, the ESS is estimated using the L2-norm of the weights:
+Given normalized (log-)weights ``w_{1:n}``, the ESS is estimated using the L2-norm of the
+weights:
 
 ```math
 \\mathrm{ESS}(w_{1:n}) = \\frac{r_{\\mathrm{eff}}}{\\sum_{i=1}^n w_i^2}
 ```
 
-where ``r_{\\mathrm{eff}}`` is the relative efficiency of the `log_weights`.
-
-    ess_is(result::PSISResult; bad_shape_nan=true)
-
-Estimate ESS for Pareto-smoothed importance sampling.
-
-!!! note
-
-    ESS estimates for Pareto shape values ``k > 0.7``, which are unreliable and misleadingly
-    high, are set to `NaN`. To avoid this, set `bad_shape_nan=false`.
+where ``r_{\\mathrm{eff}}`` is the relative efficiency of the (log-)weights.
 """
-ess_is
-
-function ess_is(r::PSISResult; bad_shape_nan::Bool=true)
-    neff = ess_is(r.weights; reff=r.reff)
-    return _apply_nan(neff, r.tail_dist; bad_shape_nan=bad_shape_nan)
-end
-function ess_is(weights; reff=1)
-    dims = _sample_dims(weights)
-    return reff ./ dropdims(sum(abs2, weights; dims=dims); dims=dims)
-end
-
-function _apply_nan(neff, dist; bad_shape_nan)
-    bad_shape_nan || return neff
-    k = pareto_shape(dist)
-    (isnan(k) || k > 0.7) && return oftype(neff, NaN)
-    return neff
+function ess_is(log_weights::AbstractArray, r_eff)
+    dims = _param_dims(log_weights)
+    return ess_is.(eachslice(log_weights; dims), r_eff)
 end
-function _apply_nan(ess::AbstractArray, tail_dist::AbstractArray; kwargs...)
-    return map(ess, tail_dist) do essᵢ, tail_distᵢ
-        return _apply_nan(essᵢ, tail_distᵢ; kwargs...)
+function ess_is(log_weights::AbstractVecOrMat, r_eff)
+    lw_max = maximum(log_weights)
+    T = typeof(exp(zero(typeof(lw_max))))
+    s1, s2 = reduce(log_weights; init=(zero(T), zero(T))) do (s1, s2), lw
+        d = lw - lw_max
+        return (s1 + exp(d), s2 + exp(2 * d))
     end
+    return r_eff * s1^2 / s2
 end
diff --git a/src/recipes/plots.jl b/src/recipes/plots.jl
deleted file mode 100644
index 04ca79dd..00000000
--- a/src/recipes/plots.jl
+++ /dev/null
@@ -1,81 +0,0 @@
-# Plots.jl recipes
-
-"""
-A module defining [`paretoshapeplot`](@ref) for plotting Pareto shape values with Plots.jl
-"""
-module PSISPlots
-
-using ..PSIS
-using RecipesBase: RecipesBase
-
-"""
-    paretoshapeplot(values; showlines=false, ...)
-    paretoshapeplot!(values; showlines=false, kwargs...)
-
-Plot shape parameters of fitted Pareto tail distributions for diagnosing convergence.
-
-`values` may be either a vector of Pareto shape parameters or a [`PSIS.PSISResult`](@ref).
-
-If `showlines==true`, horizontal lines indicating relevant Pareto shape thresholds are
-drawn. See [`PSIS.PSISResult`](@ref) for an explanation of the thresholds.
-
-All remaining `kwargs` are forwarded to the plotting function.
-
-See [`psis`](@ref), [`PSISResult`](@ref).
-
-# Examples
-
-```julia
-using PSIS, Distributions, Plots
-proposal = Normal()
-target = TDist(7)
-x = rand(proposal, 1_000, 100)
-log_ratios = logpdf.(target, x) .- logpdf.(proposal, x)
-result = psis(log_ratios)
-paretoshapeplot(result)
-```
-
-We can also plot the Pareto shape parameters directly:
-
-```julia
-paretoshapeplot(result.pareto_shape)
-```
-
-We can also use `plot` directly:
-
-```julia
-plot(result.pareto_shape; showlines=true)
-```
-"""
-paretoshapeplot, paretoshapeplot!
-
-RecipesBase.@userplot ParetoShapePlot
-
-RecipesBase.@recipe function f(plt::ParetoShapePlot; showlines=false)
-    showlines && RecipesBase.@series begin
-        seriestype := :hline
-        primary := false
-        linestyle --> [:dot :dashdot :dash :solid]
-        linealpha --> 0.7
-        linecolor --> :grey
-        y := [0 0.5 0.7 1]
-    end
-    title --> ""  # no title unless specified by the user
-    yguide --> "Pareto shape"
-    seriestype --> :scatter
-    arg = first(plt.args)
-    k = arg isa PSIS.PSISResult ? PSIS.pareto_shape(arg) : arg
-    return (PSIS.as_array(PSIS.missing_to_nan(k)),)
-end
-
-# plot PSISResult using paretoshapeplot if seriestype not specified
-RecipesBase.@recipe function f(result::PSISResult)
-    if haskey(plotattributes, :seriestype)
-        k = PSIS.as_array(PSIS.missing_to_nan(PSIS.pareto_shape(result)))
-        return (k,)
-    else
-        return ParetoShapePlot((result,))
-    end
-end
-
-end # module PSISPlots
diff --git a/src/utils.jl b/src/utils.jl
index 3b3cdad2..28029102 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -22,15 +22,15 @@ _param_axes(x::AbstractArray) = map(Base.Fix1(axes, x), _param_dims(x))
 # iterate over all parameters; combine with _selectparam
 _eachparamindex(x::AbstractArray) = CartesianIndices(_param_axes(x))
 
+function _sample_param_sizes(x::AbstractArray)
+    ndraws = size(x, 1)
+    nchains = size(x, 2)
+    nparams = prod(length, _param_axes(x); init=1)
+    return (ndraws, nchains, nparams)
+end
+
 # view of all draws for a param
 function _selectparam(x::AbstractArray, i::CartesianIndex)
     sample_dims = ntuple(_ -> Colon(), ndims(x) - length(i))
     return view(x, sample_dims..., i)
 end
-
-function _maybe_log_normalize!(x::AbstractArray, normalize::Bool)
-    if normalize
-        x .-= LogExpFunctions.logsumexp(x; dims=_sample_dims(x))
-    end
-    return x
-end
diff --git a/test/core.jl b/test/core.jl
index 77a0132b..d348d236 100644
--- a/test/core.jl
+++ b/test/core.jl
@@ -10,81 +10,44 @@ using DimensionalData: Dimensions, DimArray
 @testset "PSISResult" begin
     @testset "vector log-weights" begin
         log_weights = randn(500)
-        log_weights_norm = logsumexp(log_weights)
-        tail_length = 100
-        reff = 2.0
-        tail_dist = PSIS.GeneralizedPareto(1.0, 1.0, 0.5)
-        result = PSISResult(log_weights, reff, tail_length, tail_dist, false)
+        ess_is = 300.0
+        pareto_khat = 0.5
+        result = PSISResult(log_weights, pareto_khat, ess_is)
         @test result isa PSISResult{Float64}
-        @test issetequal(
-            propertynames(result),
-            [
-                :log_weights,
-                :nchains,
-                :ndraws,
-                :normalized,
-                :nparams,
-                :pareto_shape,
-                :reff,
-                :tail_dist,
-                :tail_length,
-                :weights,
-            ],
-        )
         @test result.log_weights == log_weights
-        @test result.weights ≈ softmax(log_weights)
-        @test result.reff == reff
-        @test result.nparams == 1
-        @test result.ndraws == 500
-        @test result.nchains == 1
-        @test result.tail_length == tail_length
-        @test result.tail_dist == tail_dist
-        @test result.pareto_shape == 0.5
-        @test result.ess ≈ ess_is(result)
+        @test result.pareto_khat == pareto_khat
+        @test result.ess_is == ess_is
 
         @testset "show" begin
             @test sprint(show, "text/plain", result) == """
                 PSISResult with 500 draws, 1 chains, and 1 parameters
-                Pareto shape (k) diagnostic values:
+                Pareto k-hat diagnostic summary:
                                     Count       Min. ESS
-                 (-Inf, 0.5]  good  1 (100.0%)  $(floor(Int, result.ess))"""
+                 (-Inf, 0.5]  good  1 (100.0%)  $(floor(Int, ess_is))"""
         end
     end
 
     @testset "array log-weights" begin
         log_weights = randn(500, 4, 3)
-        log_weights_norm = logsumexp(log_weights; dims=(1, 2))
-        log_weights .-= log_weights_norm
-        tail_length = [1600, 1601, 1602]
-        reff = [0.8, 0.9, 1.1]
-        tail_dist = [
-            PSIS.GeneralizedPareto(1.0, 1.0, 0.5),
-            PSIS.GeneralizedPareto(1.0, 1.0, 0.6),
-            PSIS.GeneralizedPareto(1.0, 1.0, 0.7),
-        ]
-        result = PSISResult(log_weights, reff, tail_length, tail_dist, true)
+        r_eff = [0.8, 0.9, 1.1]
+        ess_is = [100.0, 101.0, 102.0]
+        pareto_khats = [0.5, 0.6, 0.7]
+        result = PSISResult(log_weights, pareto_khats, ess_is)
         @test result isa PSISResult{Float64}
         @test result.log_weights == log_weights
-        @test result.weights ≈ softmax(log_weights; dims=(1, 2))
-        @test result.reff == reff
-        @test result.nparams == 3
-        @test result.ndraws == 500
-        @test result.nchains == 4
-        @test result.tail_length == tail_length
-        @test result.tail_dist == tail_dist
-        @test result.pareto_shape == [0.5, 0.6, 0.7]
-
+        @test result.pareto_khat == pareto_khats
+        @test result.ess_is == ess_is
         @testset "show" begin
             proposal = Normal()
             target = TDist(7)
             rng = StableRNG(42)
             x = rand(rng, proposal, 100, 1, 30)
             log_ratios = logpdf.(target, x) .- logpdf.(proposal, x)
-            reff = [100; ones(29)]
-            result = psis(log_ratios, reff)
+            r_eff = [100; ones(29)]
+            result = psis(log_ratios, r_eff)
             @test sprint(show, "text/plain", result) == """
                 PSISResult with 100 draws, 1 chains, and 30 parameters
-                Pareto shape (k) diagnostic values:
+                Pareto k-hat diagnostic summary:
                                        Count       Min. ESS
                  (0.5, 0.7]  okay       2 (6.7%)   99
                    (0.7, 1]  bad        2 (6.7%)   ——
@@ -94,7 +57,7 @@ using DimensionalData: Dimensions, DimArray
     end
 end
 
-@testset "psis/psis!" begin
+@testset "psis" begin
     @testset "importance sampling tests" begin
         target = Exponential(1)
         x_target = 1  # 𝔼[x] with x ~ Exponential(1)
@@ -119,30 +82,12 @@ end
                 @test r isa PSISResult
                 logw = r.log_weights
                 @test logw isa typeof(logr)
-                @test exp.(logw) == r.weights
 
-                r2 = psis(logr; normalize=false)
-                @test !(r2.log_weights ≈ r.log_weights)
-                @test r2.weights ≈ r.weights
-
-                if length(sz) > 1
-                    @test all(r.tail_length .== PSIS.tail_length(1, 400_000))
-                else
-                    @test all(r.tail_length .== PSIS.tail_length(1, 100_000))
-                end
-
-                k = r.pareto_shape
+                k = r.pareto_khat
                 @test k isa (length(sz) < 3 ? Number : AbstractVector)
-                tail_dist = r.tail_dist
-                if length(sz) < 3
-                    @test tail_dist isa PSIS.GeneralizedPareto
-                    @test tail_dist.k == k
-                else
-                    @test tail_dist isa Vector{<:PSIS.GeneralizedPareto}
-                    @test map(d -> d.k, tail_dist) == k
-                end
+                @test r.pareto_khat == k
 
-                w = r.weights
+                w = softmax(logw; dims=1:min(2, length(sz)))
                 @test all(x -> isapprox(x, k_exp; atol=0.15), k)
                 @test all(x -> isapprox(x, x_target; atol=atol), sum(x .* w; dims=dims))
                 @test all(
@@ -152,29 +97,29 @@ end
         end
     end
 
-    @testset "reff combinations" begin
-        reffs_uniform = [rand(), fill(rand()), [rand()]]
+    @testset "r_eff combinations" begin
+        r_effs_uniform = [rand(), fill(rand()), [rand()]]
         x = randn(1000)
-        for r in reffs_uniform
+        for r in r_effs_uniform
             psis(x, r)
         end
         @test_throws DimensionMismatch psis(x, rand(2))
 
         x = randn(1000, 4)
-        for r in reffs_uniform
+        for r in r_effs_uniform
             psis(x, r)
         end
         @test_throws DimensionMismatch psis(x, rand(2))
 
         x = randn(1000, 4, 2)
-        for r in reffs_uniform
+        for r in r_effs_uniform
             psis(x, r)
         end
         psis(x, rand(2))
         @test_throws DimensionMismatch psis(x, rand(3))
 
         x = randn(1000, 4, 2, 3)
-        for r in reffs_uniform
+        for r in r_effs_uniform
             psis(x, r)
         end
         psis(x, rand(2, 3))
@@ -189,17 +134,16 @@ end
                 psis(logr, rbad)
             end
             msg = String(take!(io))
-            @test occursin("All values of `reff` should be finite, but some are not.", msg)
+            @test occursin("All values of `r_eff` should be finite, but some are not.", msg)
         end
 
         io = IOBuffer()
         logr = randn(5)
         result = with_logger(SimpleLogger(io)) do
-            psis(logr; normalize=false)
+            psis(logr)
         end
         @test result.log_weights == logr
-        @test isnan(result.tail_dist.σ)
-        @test isnan(result.pareto_shape)
+        @test isnan(result.pareto_khat)
         msg = String(take!(io))
         @test occursin(
             "Warning: 1 tail draws is insufficient to fit the generalized Pareto distribution.",
@@ -215,11 +159,10 @@ end
             vcat(ones(50), fill(-Inf, 435)),
         ]
             result = with_logger(SimpleLogger(io)) do
-                psis(logr; normalize=false)
+                psis(logr)
             end
             @test skipnan(result.log_weights) == skipnan(logr)
-            @test isnan(result.tail_dist.σ)
-            @test isnan(result.pareto_shape)
+            @test isnan(result.pareto_khat)
             msg = String(take!(io))
             @test occursin("Warning: Tail contains non-finite values.", msg)
         end
@@ -229,58 +172,48 @@ end
         x = rand(rng, Exponential(50), 1_000)
         logr = logpdf.(Exponential(1), x) .- logpdf.(Exponential(50), x)
         result = with_logger(SimpleLogger(io)) do
-            psis(logr; normalize=false)
+            psis(logr)
         end
         @test result.log_weights != logr
-        @test result.pareto_shape > 0.7
+        @test result.pareto_khat > 0.7
         msg = String(take!(io))
-        @test occursin(
-            "Warning: Pareto shape k = 0.72 > 0.7. $(PSIS.BAD_SHAPE_SUMMARY)", msg
-        )
+        @test occursin("Warning: Pareto k-hat = 0.72 > 0.7. $(PSIS.BAD_KHAT_SUMMARY)", msg)
 
         io = IOBuffer()
         with_logger(SimpleLogger(io)) do
-            PSIS.check_pareto_shape(PSIS.GeneralizedPareto(0.0, 1.0, 1.1))
+            PSIS.check_pareto_khat(1.1)
         end
         msg = String(take!(io))
         @test occursin(
-            "Warning: Pareto shape k = 1.1 > 1. $(PSIS.VERY_BAD_SHAPE_SUMMARY)", msg
+            "Warning: Pareto k-hat = 1.1 > 1. $(PSIS.VERY_BAD_KHAT_SUMMARY)", msg
         )
 
         io = IOBuffer()
         with_logger(SimpleLogger(io)) do
-            PSIS.check_pareto_shape(PSIS.GeneralizedPareto(0.0, 1.0, 0.8))
+            PSIS.check_pareto_khat(0.8)
         end
         msg = String(take!(io))
-        @test occursin(
-            "Warning: Pareto shape k = 0.8 > 0.7. $(PSIS.BAD_SHAPE_SUMMARY)", msg
-        )
+        @test occursin("Warning: Pareto k-hat = 0.8 > 0.7. $(PSIS.BAD_KHAT_SUMMARY)", msg)
 
         io = IOBuffer()
         with_logger(SimpleLogger(io)) do
-            PSIS.check_pareto_shape(PSIS.GeneralizedPareto(0.0, 1.0, 0.69))
+            PSIS.check_pareto_khat(0.69)
         end
         msg = String(take!(io))
         @test isempty(msg)
 
-        tail_dist = [
-            PSIS.GeneralizedPareto(0, NaN, NaN),
-            PSIS.GeneralizedPareto(0, 1, 0.69),
-            PSIS.GeneralizedPareto(0, 1, 0.71),
-            PSIS.GeneralizedPareto(0, 1, 1.1),
-        ]
+        khats = [NaN, 0.69, 0.71, 1.1]
         io = IOBuffer()
         with_logger(SimpleLogger(io)) do
-            PSIS.check_pareto_shape(tail_dist)
+            PSIS.check_pareto_khat(khats)
         end
         msg = String(take!(io))
         @test occursin(
-            "Warning: 1 parameters had Pareto shape values 0.7 < k ≤ 1. $(PSIS.BAD_SHAPE_SUMMARY)",
+            "Warning: 1 parameters had Pareto k-hat values in (0.7, 1]. $(PSIS.BAD_KHAT_SUMMARY)",
             msg,
         )
         @test occursin(
-            "Warning: 1 parameters had Pareto shape values k > 1. $(PSIS.VERY_BAD_SHAPE_SUMMARY)",
-            msg,
+            "Warning: 1 parameters had Pareto k-hat > 1. $(PSIS.VERY_BAD_KHAT_SUMMARY)", msg
         )
         @test occursin(
             "Warning: For 1 parameters, the generalized Pareto distribution could not be fit to the tail draws.",
@@ -298,17 +231,22 @@ end
         logr = permutedims(logr, (2, 3, 1))
         @testset for r_eff in (0.7, 1.2)
             r_effs = fill(r_eff, sz[1])
-            result = @inferred psis(logr, r_effs; normalize=false)
+            result = @inferred psis(logr, r_effs)
             logw = result.log_weights
             @test !isapprox(logw, logr)
-            basename = "normal_to_cauchy_reff_$(r_eff)"
+            basename = "normal_to_cauchy_r_eff_$(r_eff)"
             @test_reference(
                 "references/$basename.jld2",
-                Dict("log_weights" => logw, "pareto_shape" => result.pareto_shape),
+                Dict(
+                    "log_weights" => logw,
+                    "pareto_khat" => result.pareto_khat,
+                    "ess_is" => result.ess_is,
+                ),
                 by =
                     (ref, x) ->
                         isapprox(ref["log_weights"], x["log_weights"]; rtol=1e-6) &&
-                            isapprox(ref["pareto_shape"], x["pareto_shape"]; rtol=1e-6),
+                        isapprox(ref["pareto_khat"], x["pareto_khat"]; rtol=1e-6) &&
+                        isapprox(ref["ess_is"], x["ess_is"]; rtol=1e-6)
             )
         end
     end
@@ -336,8 +274,8 @@ end
             result = @inferred psis(logr)
             @test result.log_weights isa DimArray
             @test Dimensions.dims(result.log_weights) == Dimensions.dims(logr)
-            for k in (:pareto_shape, :tail_length, :tail_dist, :reff)
-                prop = getproperty(result, k)
+            @testset for k in (:pareto_khat, :ess_is)
+                prop = getfield(result, k)
                 @test prop isa DimArray
                 @test Dimensions.dims(prop) == Dimensions.dims(logr, (:param,))
             end
diff --git a/test/ess.jl b/test/ess.jl
index 05c7fcd3..58ad17cd 100644
--- a/test/ess.jl
+++ b/test/ess.jl
@@ -4,39 +4,21 @@ using PSIS
 using Test
 
 @testset "effective sample size" begin
-    w = fill(0.01, 100)
-    reff = rand()
-    @test ess_is(w) ≈ 100
-    @test ess_is(w; reff=reff) ≈ 100 .* reff
+    logw = fill(randn(), 100)
+    @test PSIS.ess_is(logw, 1) ≈ 100
+    r_eff = rand()
+    @test PSIS.ess_is(logw, r_eff) ≈ 100 * r_eff
 
-    w = zeros(100)
+    logw = fill(-Inf, 100)
     i = rand(1:100)
-    w[i] = 1
-    @test ess_is(w) ≈ 1
-    @test ess_is(w; reff=reff) ≈ 1 .* reff
-
-    logw = randn(100)
-    result = PSISResult(logw, 1.5, 20, PSIS.GeneralizedPareto(0.0, 1.0, 0.6), false)
-    @test ess_is(result) ≈ ess_is(result.weights; reff=1.5)
-
-    result = PSISResult(logw, 1.5, 20, PSIS.GeneralizedPareto(0.0, 1.0, 0.71), false)
-    @test isnan(ess_is(result))
-    @test ess_is(result; bad_shape_nan=false) ≈ ess_is(result.weights; reff=1.5)
+    logw[i] = 0
+    @test PSIS.ess_is(logw, 1) ≈ 1
+    @test PSIS.ess_is(logw, r_eff) ≈ r_eff
 
     logw = randn(100, 4, 3)
-    tail_dists = [
-        PSIS.GeneralizedPareto(0.0, 1.0, 0.69),
-        PSIS.GeneralizedPareto(0.0, 1.0, 0.71),
-        PSIS.GeneralizedPareto(0.0, NaN, NaN),
-    ]
-    reff = [1.5, 0.8, 1.0]
-    result = PSISResult(logw, reff, [20, 20, 20], tail_dists, false)
-    ess = ess_is(result)
-    @test ess isa Vector
-    @test length(ess) == 3
-    @test ess[1] ≈ ess_is(result.weights; reff=reff)[1]
-    @test isnan(ess[2])
-    @test isnan(ess[3])
-    ess = ess_is(result; bad_shape_nan=false)
-    @test ess ≈ ess_is(result.weights; reff=reff)[1:3]
+    w = softmax(logw; dims=(1, 2))
+    dims = (1, 2)
+    @test PSIS.ess_is(logw, r_eff) ≈ r_eff ./ dropdims(sum(abs2, w; dims); dims)
+    r_eff = rand(3)
+    @test PSIS.ess_is(logw, r_eff) ≈ r_eff ./ dropdims(sum(abs2, w; dims); dims)
 end
diff --git a/test/plots.jl b/test/plots.jl
deleted file mode 100644
index 5f437edb..00000000
--- a/test/plots.jl
+++ /dev/null
@@ -1,59 +0,0 @@
-using PSIS
-using Plots
-using Test
-
-@testset "PSISPlots" begin
-    @testset "$f" for f in (PSISPlots.paretoshapeplot, PSISPlots.paretoshapeplot!),
-        sz in (100, (100, 10))
-
-        result = psis(randn(sz...))
-
-        @testset for values in (result, result.pareto_shape)
-            plot()
-            plt = f(values)
-            @test plt isa Plots.Plot
-            @test length(plt.series_list) == 1
-            @test plt[1][1][:x] == Base.OneTo(result.nparams)
-            @test plt[1][1][:y] == PSIS.as_array(result.pareto_shape)
-            @test plt[1][1][:seriestype] == :scatter
-            @test plt[1][:yaxis][:guide] == "Pareto shape"
-        end
-
-        plot()
-        plt2 = f(result; showlines=true)
-        @test plt2 isa Plots.Plot
-        @test length(plt2.series_list) == 5
-        linestyles = [:dot, :dashdot, :dash, :solid]
-        yvals = [0, 0.5, 0.7, 1]
-        for i in 1:4
-            @test plt2[1][i][:y] == fill(yvals[i], 3)
-            @test plt2[1][i][:linestyle] == linestyles[i]
-            @test plt2[1][i][:linecolor] ==
-                RGBA(0.5019607843137255, 0.5019607843137255, 0.5019607843137255, 1.0)
-            @test plt2[1][i][:linealpha] == 0.7
-        end
-    end
-
-    @testset "plot(::PSISResult)" begin
-        result = psis(randn(100, 10))
-        for showlines in (true, false)
-            plt = PSISPlots.paretoshapeplot(result; showlines=showlines)
-            plt2 = plot(result; showlines=showlines)
-            @test length(plt.series_list) == length(plt2.series_list)
-            for (s1, s2) in zip(plt.series_list, plt2.series_list)
-                @test s1[:seriestype] === s2[:seriestype]
-                @test isequal(s1[:x], s2[:x])
-                @test isequal(s1[:y], s2[:y])
-            end
-        end
-    end
-
-    @testset "plot(::PSISResult; seriestype=:path)" begin
-        result = psis(randn(100, 2, 10))
-        plt = plot(result; seriestype=:path)
-        @test length(plt.series_list) == 1
-        @test plt[1][1][:x] == eachindex(result.pareto_shape)
-        @test plt[1][1][:y] == result.pareto_shape
-        @test plt[1][1][:seriestype] == :path
-    end
-end
diff --git a/test/references/normal_to_cauchy_reff_0.7.jld2 b/test/references/normal_to_cauchy_r_eff_0.7.jld2
similarity index 75%
rename from test/references/normal_to_cauchy_reff_0.7.jld2
rename to test/references/normal_to_cauchy_r_eff_0.7.jld2
index 1023ebad..6602d2d7 100644
Binary files a/test/references/normal_to_cauchy_reff_0.7.jld2 and b/test/references/normal_to_cauchy_r_eff_0.7.jld2 differ
diff --git a/test/references/normal_to_cauchy_reff_1.2.jld2 b/test/references/normal_to_cauchy_r_eff_1.2.jld2
similarity index 97%
rename from test/references/normal_to_cauchy_reff_1.2.jld2
rename to test/references/normal_to_cauchy_r_eff_1.2.jld2
index 5ae886ce..926f3fa5 100644
Binary files a/test/references/normal_to_cauchy_reff_1.2.jld2 and b/test/references/normal_to_cauchy_r_eff_1.2.jld2 differ
diff --git a/test/runtests.jl b/test/runtests.jl
index 47b31094..e24f320f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,5 +6,4 @@ using Test
     include("generalized_pareto.jl")
     include("core.jl")
     include("ess.jl")
-    include("plots.jl")
 end