Don't error for very low log-weights (#28)

* Don't rescale distribution

* Add basic test

* Increment version number

Project.toml

@@ -1,7 +1,7 @@
 name = "PSIS"
 uuid = "ce719bf2-d5d0-4fb9-925d-10a81b42ad04"
 authors = ["Seth Axen <[email protected]> and contributors"]
-version = "0.2.6"
+version = "0.3.0"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

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

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

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