diff --git a/Project.toml b/Project.toml index 6a424db1..128dabfb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PSIS" uuid = "ce719bf2-d5d0-4fb9-925d-10a81b42ad04" authors = ["Seth Axen and contributors"] -version = "0.2.6" +version = "0.3.0" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/src/core.jl b/src/core.jl index 7546aa6e..a64bcb87 100644 --- a/src/core.jl +++ b/src/core.jl @@ -30,7 +30,9 @@ Result of Pareto-smoothed importance sampling (PSIS) using [`psis`](@ref). - `log_weights_norm`: the logarithm of the normalization constant of `log_weights` - `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` + `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`. # Diagnostic @@ -303,7 +305,5 @@ function psis_tail!(logw, logμ, M=length(logw), improved=false) logw[i] = min(log(_quantile(tail_dist_adjusted, p[i])), 0) + logw_max end end - # undo scaling for the tail distribution - tail_dist = scale(tail_dist_adjusted, exp(logw_max)) - return logw, tail_dist + return logw, tail_dist_adjusted end diff --git a/src/generalized_pareto.jl b/src/generalized_pareto.jl index f7109e27..d2ef1f8e 100644 --- a/src/generalized_pareto.jl +++ b/src/generalized_pareto.jl @@ -184,7 +184,3 @@ function prior_adjust_shape(d::Distributions.GeneralizedPareto, n, ξ_prior=1//2 ξ = (n * d.ξ + nobs * ξ_prior) / (n + nobs) return Distributions.GeneralizedPareto(d.μ, d.σ, ξ) end - -function scale(d::Distributions.GeneralizedPareto, s) - return Distributions.GeneralizedPareto(d.μ * s, d.σ * s, d.ξ) -end diff --git a/test/core.jl b/test/core.jl index 870bed05..f647996e 100644 --- a/test/core.jl +++ b/test/core.jl @@ -264,6 +264,11 @@ end end end + # https://github.com/arviz-devs/PSIS.jl/issues/27 + @testset "no failure for very low log-weights" begin + psis(rand(1000) .- 1500) + end + @testset "compatibility with arrays with named axes/dims" begin param_names = [Symbol("x[$i]") for i in 1:10] iter_names = 101:200